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ABSTRACT 

There exists strong circumstantial evidence from their eccentric orbits that most of the known extra-solar planetary systems are the 
survivors of violent dynamical instabilities. Here we explore the effect of giant planet instabilities on the formation and survival of 
terrestrial planets. We numerically simulate the evolution of planetary systems around Sun-like stars that include three components: 
(i) an inner disk of planetesimals and planetaiy embryos, (ii) three giant planets at Jupiter-Saturn distances, and (iii) an outer disk of 
planetesimals comparable to estimates of the primitive Kuiper belt. We calculate the dust production and spectral energy distribution 
of each system by assuming that each planetesimal particle represents an ensemble of smaller bodies in coUisional equilibrium. Our 
main result is a strong correlation between the evolution of the inner and outer parts of planetary systems, i.e. between the presence of 
terrestrial planets and debris disks. Strong giant planet instabilities - that produce very eccentric surviving planets - destroy all rocky 
material in the system, including fully-formed terrestrial planets if the instabilities occur late, and also destroy the icy planetesimal 
population. Stable or weakly unstable systems allow terrestrial planets to accrete in their inner regions and significant dust to be 
produced in their outer regions, detectable at mid-infrared wavelengths as debris disks. Stars older than ~ 100 Myr with bright 
cold dust emission (in particular al A ~ 70yum) signpost dynamically calm environments that were conducive to efficient terrestrial 
accretion. Such emission is present around ~16% of billion-year old Solar-type stars. 

Our simulations yield numerous secondary results: 1) The typical eccentricities of as-yet undetected terrestrial planets are ~0.1 but 
there exists a novel class of terrestrial planet system whose single planet undergoes large amplitude oscillations in orbital eccentricity 
and inclination; 2) By scaling our systems to match the observed semimajor axis distribution of giant exoplanets, we predict that 
terrestrial exoplanets in the same systems should be a few times more abundant at ~ 0.5 AU than giant or terrestrial exoplanets at 1 
AU; 3) The Solar System appears to be unusual in teims of its combination of a rich terrestrial planet system and a low dust content. 
This may be explained by the weak, outward-directed instability that is thought to have caused the late heavy bombardment. 

Key words, planetary systems: formation — methods: n-body simulations — circumstellar matter — infrared stars — Kuiper belt 
— Solar System 
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.Cd. 1. Introduction lO'Brien et al.l2006l) . From roughly a few to a few tens of AU. gi- 

ant planet core s grow and accrete gaseous envelopes if t he condi 



ant planet core s grow ana accrete gaseous envelopes it m e conai- 
tions are right (iPollack et al.lll996l: IXiibert et al.ll2005l) . Despite 



Circumstellar disks of gas and dust are expected to pro- 
duce t hree broad classes o f planets in radially-segregated their lai-ge masses, gas giants must form wit hin the few million 
zones (iKokubo & Idal l2()03) . The inner disk forms terres- vear lifetime of gaseous disks (iHaisch et al. 2001 : Pascucci et aU 
trial (rocky) planets because it contains too litde solid mass |2006|; [Kennedy & Kenyon|l200g) and be present during the late 
to rapidly accrete giant planet cores, which are thought to phases of terrestrial planet growth. Resonant interactions, aris- 
form preferentially beyond the snow line where the sur- ingnotjust from giant planets but also from the changing surface 
face density in solids is higher b ecause of ice c ondensation density of the gas disk itself (Nagasawa et al. 2005), thus likely 
and the isolation mass is larger (iLissaued [T9871) . TeiTestrial play a role in teiTestrial planet formation. Finally, in the outer 
planets form in 10-100 million years (Myr) via collisional regions of planetai-y systems the growth time scale exceeds the 
agglomeration of Moon- to Mars-sized planetary embryos lifetime of the gas disk, and the end point of accretion is a belt 
and a swarm of 1-10^ km s i zed planetesima ls ( Chambers of Pluto-sized (and smaller) bodies (jKenyon & Luu||1998|). 
1200 It iKenvon & Bromlevll200d: iRavmond et alJl2006a i2009a 
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Our understanding of planetary formation in these different 
zones is constrained by a variety of observations. The initial 
conditions of inner disks can be probed by observations of hot 
dust. Such observations have show n evidence for grain gr owth 
in individual protoplanetary disks (iBouwman et all 12008 1). and 
larger samples have shown that hot dust disappears from these 
disks within a few Myr and that cooler dust takes lo nger to 
disap pear (Haisch et al. 2001; Mamajek et al. 2004; Me ver et alj 
l2008h . Hundreds of close-in planets - the outcome of planet for- 
mation - have been detected by radial velocity and transit obser- 
vations, with masses down to a few (e.g., [B utler et al. 2006; 
lUdrv et al. 120071; iLeger et al .'2009: Batalha et al. 2011). The fre- 
quency of planets on short-period orbits has been found to be a 
function of the planet mass; less massive planets are signi ficantly 
more com mon than high-mass "hot Jupiters" ( Mavor et al.l2009l: 
iHoward et al. 2010). The frequency of 3-10 planets has been 
measured at 12%, and by extrapolation the frequ ency of 0.5-2 
Mffi i s 23% (for orbital periods less than 50 davs; iHoward et all 
I2010h . 

Despite the existence of hot Jupiters around ~ 1 % of solar- 
type stars (Howard et al. 2010), the vast majority of g iant plan- 
ets i s found beyond 1 AU (iButler et al.i 120061; iUdrv & SantosI 
|2007|) . The absolute frequency of gi ant planets is poorly c on- 
strained; estimate s range from 10% (ICumming et al.ll2008l) to 
more than 50% (iGould et al.l |2.010'). These planets are char- 
acterized by their broad eccentricity distribu tion that includes 
several planets with e > 0.9 But ler et al.l ([2006). This dis- 
tribution is quantitatively reproduced if giant planets form in 
systems with multiple planets, and dynamical instabilities oc- 
curr ed in 70-100 % of all observed systems (IChatteriee et al.l 
[2008; June & Tremainel l2008 : Raymond.etalJ[2010|). The onset 
of instability may be caused by the chang ing planet-planet sta - 
bility criterion as the gas disk dissipates ("iwa saki et al.ll200"l1) . 
resonant migration (Adams & Laughlin 2003), or chaotic dy- 
namics ( Chambers et al.lll99 6). Whatever the trigger, the insta- 
bihty leads to a phase of planet-planet scattering and in most 
cases to the eventual removal of one or more planets from 
the system by collision or hyperbolic ejection (iRasio & FordI 
|l996; Weidenschilling & Marzari 1996). It is the surviving plan- 
ets that match the observed distribution. The properties of the 
outer Solar System are consistent with the gi ant planets hav- 
ing f ormed in a similarly unstable configuration (iThommes et alJ 
119991) . though the low eccentricities of Jupiter and Saturn 
(ej «; eg s; 0.05), and the lat e timing of the Late Heavy 
Bombardment dStrom et al.ll2005h . hint that the instability that 
occurred in our Solar Syst em was weak and occu rred too late to 
affect terrestrial accretion (iMorbidelli et al.|[20Tol ). 

Sub-millimeter observations of dust disks around young 
stars provide information on the initial conditions in outer 
planet-forming disks, and by extrapolation to entire disks. These 
observations suggest that the typi cal protoplan etary disk has 
a mass of 0.001-0.1 Solar masses ([Andrews & Williams 2005, 
l2007at lEisne r et"ar 2008) and a radial s urface density profile 
of roughly r'C' s-') (MiindvetalJ 120001; [Andrews & Williams 
I2007bt1lsella"etan 120101) . There appears to be a roughly Un- 
ear trend between the dust mass and the stellar mass such that 
the typical p rotoplanetary disk contains a few percent of the 
stellar mass ( Andrews & Williamsll2007a[) . which has important 
impli cations for the planet frequency as a function of st ellar 
mass (llda & Lid[2005l;lRavmond"etan[2007t [G reave 

Debris disks - warm or cold dust observed around older 
stars, typically at infrared wavelengths (/I ~ 10 -lOOyum)- probe 
outer disks after planet formation has completed. Debris disks 
provide evidence for the existence of leftover planetesimals be- 



cause the lifetime of observed dust particles under the effects 
of collisions and radiation forces is far shorter than the typi- 
cal stellar age, implying a replenishment via collisional grind- 
ing of larger bodies. Spitzer observations show that about 15% 
of solar-type stars younger than 300 Myr have significant dust 
excesses at 24fj.m but that t his fraction decrea s es to about 3% 
for stars older than 1 Gyr (Mever et al."2008', ICarpenter et al.l 
2009). At 70yum the fraction of stars showing significant excesses 
is roughly constant in time (at ~ 16%) but the u pper envelope 
of the di stribution decrea s es with the stellar age (iTrilling et al.l 
2008; C aspar et al.l 120091; ICarpenter et al1l2009l) . A stars have 
a significantly higher fraction of dust excesses at both 24 and 
70jum and the excesses themselves a re brighter tha n for FGK 
stars but the dust lifetime is shorter (iSu et al.ll2006l) . Very few 
debris disks are known around M dwarfs and the d ust brightness 
is nec essarily far lower than for higher-mass stars (ICautier et al.l 
[2OOI . 

Here we perform a large ensemble of N-body simulations to 
model the interactions between the different radial components 
of planetary systems: formation and survival of terrestrial plan- 
ets, dynamical evolution and scattering of giant planets, and dust 
production from collisional grinding. By matching the orbital 
distribution of the giant planets, we infer the characteristics of 
as-yet undetected terrestrial planets in those same systems. We 
post-process the simulations to compute the spectral energy dis- 
tribution of dust by treating planetesimal particles as aggregates 
in collisional equilibrium (Dohnanvi 1969) to calculate the inci- 
dent and re-emitted flux (.Booth et al.. 2009) . We then correlate 
outcomes in the different radial zones and link to two key ob- 
servational quantities: the orbital properties of giant planets and 
debris disks. 

The paper is structured as follows. In section 2 we explain 
our choice of initial conditions, the integration method, and our 
debris disk calculations. In section 3 we demonstrate the detailed 
evolution of a single simulation in terms of its dynamics, the for- 
mation of terrestrial planets, and the dust evolution. In section 4 
we present results from our fiducial set of 152 simulations and 
explore the correlations between terrestrial planet formation ef- 
ficiency, giant planet characteristics and debris disks. In Section 
5 we discuss the implications of our results: we scale our simu- 
lations to the known eccentricity-semimajor axis distribution of 
giant planets, discuss the observed debris disks in known exo- 
planet systems, and evaluate the Solar System in the context of 
our results. We conclude in Section 6. 

In a companion paper (Raymond et al 201 1 ; referred to in the 
text as Paper 2) we test the effects of several system parameters 
on this process: the giant planet masses and mass distribution, 
the width of the outer planetesimal disk, the mass distribution of 
the outer planetesimal disk, and the presence of disk gas at the 
time of giant planet instabilities. We also calibrate our results to 
the statistics of known debris disks to produce an estimate of the 
fraction of solar-type stars that host terrestrial planets {rjEarth in 
the Drake equation). 

2. Methods 

We use numerical simulations to study the formation of terres- 
trial planets from massive embryos, the dynamical evolution of 
fully-formed gas and ice giants, and the long-term evolution of 
debris disks, starting from an assumed set of initial conditions 
that are specified at the epoch when the protoplanetary gas disk 
dissipates. The goal is to study how these processes, occurring 
in spatially distinct regions of the planetary system, are cou- 
pled, and to quantify the expected dispersion in the final out- 
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come that arises from the chaotic nature of the evolution. In what 
follows, we describe in detail the specific initial conditions that 
we adopt, together with the numerical methods. It is worth em- 
phasizing at the outset that our initial conditions are not unique. 
Observationally, little is known about the radial distribution of 
planetesimal formation within protoplanetary disks, and as a re- 
sult the single most important initial condition for subsequent 
planet formation is not empirically well-constrained. Our fun- 
damental assumption is that planetesimals form across a wide 
range of disk radii (extending out, in particular, to encompass 
a cold outer debris disk), with a smooth radial distribution. We 
further assume that cores are able to form quickly enough to 
yield planets at least as massive as ice giants in typical disk^. 
Given these assumptions, it is likely that the system-to-system 
variation in the masses of giant planets - arising ultimately from 
dispersion in the masses and radii of the protoplanetary disks - 
exceeds that of the terrestrial planets or that of the debris disk, 
because small variations in the formation time scale of cores re- 
sult in large changes to the final envelope mass. 

Our models are an extension of what might be described as 
the "classical" scenario for planet formation. Substantially dif- 
ferent scenarios are also possible. In particular , planetesimal for- 
mation is poorly understood (for a review, see thian^ & YoudinI 
[2010), and may be coupled to the level of intrinsic turbulence 
within the gaseous disk, whic h probably varies with radius 
(lGammielll996tlArmitagell201 ih . If planetesimal formation effi- 
ciency varies dramatically with radius, the initial conditions for 
subsequent planet formation would be radically different from 
what we have assumed, leading to qualitatively different con- 
clusions. For example, it could be true that the zones of terres- 
trial and giant planet formation are typically dynamically well- 
separated, contrary to what is implied by our initial conditions. 
In such a model, the coupling between the giant and terrestrial 
planets would be much weaker than occurs given our assump- 
tions, such that only exceptionally violent giant planet scattering 
- such as occurs when instabilities drive e ^ 1 - would affect 
the terrestrial zone. Similar caveats apply to the outer debris disk 
region. 



2. 1 . Initial conditions 

The initial conditions for planet formation are expected to vary 
with stellar mass, due both to the relatively well-understood 
variation of the location of the snow line with stellar type 
(iKennedv & Kenvorj [200i) . and as a r esult of any systematic 
variat ions in the star-to-disk mass ratio (lAlexander & Armitagd 
l2006l) . We consider Solar-mass stars, and start our simulations 
from highly idealized initial conditions that represent the pre- 
dicted state of a planetary system at the time of the dissipation of 
the gaseous protoplanetary disk. There are three radial zones: an 
inner zone containing 9M^ in planetary embryos and planetesi- 
mals from 0.5 to 4 AU, three giant planets on marginally stable 
orbits from Jupiter's orbital distance of 5.2 AU out to ~10 AU 
(depending on the masses), and an outer 10 AU-wide disk of 
planetesimals containing 50 M®. In the majority of simulations 
the giant planets underwent a violent phase of scattering but in a 
significant fraction they did not. Figure [T] shows an example set 
of initial conditions. 
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' Recent calculati ons of planet formatio n via core accretion support 
such an assumption jMovshovitz et al.l201 0). though it should be noted 
that these are subject to substantial uncertainties bec ause the physics o f 
Type I migration remains imperfectly understood (iLubow & Idal2010h . 



Fig. 1. An example of the initial conditions assumed to exist 
shortly after the epoch of gas disk dispersal. The terrestrial 
planet formation zone contains planetary embryos and lower 
mass planetesimals, while the outer planetesimal belt contains 
a population of relatively low-mass (compared to the giant plan- 
ets) bodies that are represented numerically by equal mass 
particles. The intermediate zone contains 3 fully-formed giant 
planets whose spacing is such that they are marginally unstable 
against dynamical instability. 



2.1.1. Terrestrial planet-forming region 

The terrestrial planet forming region extends from 0.5 to 4 AU. 
As in the bulk of terrestrial planet formation simulations the 
inner boundary was chosen as a compromise between simula- 
tion run time and capturing the dynamics of the inner disk (e.g. 
ICham bers'2001). The outer boundary corresponds to the approx- 
imate stability boundary with the innermost giant planet (at 5.2 
AU), depending on that planet's mass, such that obj ects beyond 
this limit woul d be im mediately destabilized ( Marchal & BozisI 
1982; Gladmanl ll993h . Within this zone, we adopt initial condi- 
tions that are motivate d by "chao tic growth" models, of the type 
su mmarized b y Goldre ich et al.l ( 1200 4) and calculated in detail 
bv lKenvon & Bromlev ( l2006l) . Specificallv. the terrestrial zone 
contains a total of 9 in 49 planetary embryos and 500 plan- 
etesimals, with equal mass in each component. The embryo mass 
of 0.09 Me is a factor of a few to ten larger than t hat used in 
recent simulations of late-st age terrestrial accretion (ICharnbersI 
boOlLlRavmond et al. 2006a. l2009d:lMorishima et al.l2010l) . but 
is within a factor of ~2 of that obtained from semi-analytic 
models of oligarchic growth (Chambers 2006) over the radial 
range between 1 AU and 3 AU at 3 Myr Of course, the out- 
come of late-stage accretion has been shown to be insensitive 
to the embryo mass distribution but rather to d epend mainly 
on th e large-scale mass distribution of the disk (iKokubo et al.l 
l 2006h . The mass is distributed following a radial surface density 
profile S oc r ^ This profile is consistent with an inward ex - 
trapolation of that measured for DG Tau by llsella et aP (1201 Ol) . 
who infer a slope for the disk surface density via interferometric 
observations of thermal dust emission at mm wavelengths. We 
do not consider the possibility of variations in the disk mass, 
due either to dispersion in the surface density of the proto- 
planetary disk at early times, or arising from dynamical effects 
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such as the passag e of massive planets migrating inward due to 
Type II migration ('Fogg & Nelsonll2007LlRavmond et al.ll2006 a: 
[Mandell et al. 2007). Nor do we consider the effect of the chang- 
ing gravitational potential of the dissipating ga s disk, which in 
some cases could influe nce terrestrial accretion (iNagasawa et al.l 
120051: iMorishima et all 2010). Modeling these effects would re- 
quire starting our calculations at an early epoch, when the gas 
disk was still present, and would introduce additional physical 
uncertainties. Our disk is comparable in mass to the "minimum - 
mass solar nebula" (MMSN) model dWeidenschiUind [TqTTI) . 
with 2.6 [5.1] Me inside 1.5 [2.5] AU, though we do not adopt 
the MMSN surface density slope. 

The embryo mass increases only slightly with orbital dis- 
tance , in agreement with models for embryo growth (IChambersI 
l2006 l). In practice, embryos are spaced by A mutual Hill radii 

RH,m- 

1 .lmi+m2\^'^ ... 

I^Hm — —yci\ + a2)\ , (1) 

• 2^ '\ / ^ ' 

where a\ and 02 are orbital distances of two adjacent embryos, 
m\ and 1112 are their masses, and is the stellar mass. We chose 
A to lie between 8-16 but decreasing systematically with orbital 
distance to avoid large variations in embryo mass in different re- 
gions of the disk. In practice, we used A = 8 -1- 61a, where 5 
is randomly chosen between zero and 8 and a is the orbital dis- 
tance in AU. The mass resolution of the terrestrial component 
of our model is within a factor of 2-4 of the best current sim- 
ulations of terrestrial p lanet formation (iRavmond et al.ll2009ct 
IMorishima et al.ll2010l) . Embryos were given randomly-chosen 
initial eccentricities of up to 0.02 and initial inclinations of up to 
0.1°. Planetesimals were given initial eccentricities of up to 0.02 
and initial inclinations of up to 1°. 

2.1.2. Giant Planets 

Our giant planet distribution is motivated by observations of 
extrasolar planetary systems, which show that the giant planet 
population detecte d at a < 5 A U exhibits a broad eccentricity 
distribution (Butl er et al.ll2006l) . The ecc entricity distibution i s 
not strongly affected by selection effects (Gumming et al.ll2008h . 
which do however dominate the radial distribution (which is 
essentially unknown beyond about 5 AU). There is no unique 
theoretical interpretation of this result, but the most common 
explanation is that the eccentricity represents the endpoint of 
a dy namical scattering phase. Models built u pon this assump- 
tion (lJuric & Tremaine 2008; C hatteriee et al.l l2008) can repro- 
duce quantitatively the observed f{e) for extrasolar planets. 
When the effects of outer planetesimal disks are included, they 
are also consistent with near-circular orbits bein g typical for low 
mass giant planet systems at larger orbital radii (iRavmond et al.l 
[2010, as is found in the Solar System). 

Here, we assume that marginally unstable initial condi- 
tions are also the rule for orbital separations modestly larger 
than those needed to explain the current observations of mas- 
sive extrasolar planets. This may not be true, but it repre- 
sents the simplest extrapolation of current observational re- 
sults. We model systems of three giant planets, with the inner- 
most giant planet given an orbital semimajor axis of 5.2 AU, 
the same as Jupiter's. This is an arbitrary choice which al- 
lows for easy comparison with the Solar System (assuming 
that Jupiter formed at 5.2 AU). Two additional planets are 
spaced outward with randomly-chosen separations of (in the 
fiducial runs) 4 - SRn.m- This separation was chosen to select 



for systems that will likely become unstable on a timescale 
of 100,000 years or longer (iMarzari &Weidenschillingl 120021 
Chatteriee et al. 2008). This value of 10^ years is the expected 
timescale for the final dissipation of the gaseous p rotoplane- 
tary disk dWolk & Waltedll996t lErcolano et all|201 ll) . although 
substantial structural change s to the disk occur over longer 
timescales (C urrie et al.l2009l) . Because damping within the disk 
tends to stabilize planetary systems Jlwasaki et a l. 2002), we ex- 
pect the natural instability timescale to be on the order of the gas 
dissipation timescale. 

The outcome of scattering is dependent on both the ini- 
tial mass distribution (whic h is related to the well-constrained 
observed mass distribution iButler et al 1 l2006h . and on the dis- 
tributi on of mass ratios between planets in individual sys- 
tems (iFord et al.1 120031: IRavmond et al] \2Dl(f). For the results 
shown in the paper we focus on a fiducial set of runs (called 
mixed in paper 2), in which the planet masses are chosen to fol- 
low the observed distribution of extra-solar planets. 



where masses are chosen between one Saturn mass and 3 Jupiter 
masses. This range is chosen so as to include the majority of 
well-observed exoplanet systems (more massive planets are rel- 
atively uncommon, while selection effects become stronger for 
very low mass giants), but there is nothing particularly special 
about our choice. Similar models, that fit the observed con- 
straints equally well, could almost certainly be constructed as- 
suming different mass ranges, although it might be necessary 
to adopt different numb ers of planets or to as sume correlations 
between their masses (iRavmond et al.l[2010l) . In our runs, the 
masses of individual planets are chosen independently. 

The radii of giant planets affect their early dynamics, by al- 
tering the ratio of physical collisions to scattering events. We 
adopt a constant bulk density, p = 13 g cm"^, independent of 
mass, that matches that of the (present-day) Jupiter. Since young 
giant planets will assuredly have larger radii, our assumption un- 
derestimates the true number of physical collisions. We do not, 
however, consider this to be a major source of error Planets in 
our mass range, formed via core accretion, are only modestly 
inflated at ages as short as a Myr (iMarlev et al.ll200'7b . and col- 
lisions are suppressed at a > 5 AU because of the relatively 
modest orbital velocities at these radii (Ford et al. 2001). 

As noted above, the initial conditions that we adopt for 
the giant planets are motivated by strictly empirical conditions: 
the observed mass function of extrasolar planets, and the re- 
quirement (within the context of a scattering model) that most 
multiple planet systems are eventually dynamically unstable. 
However, they are also broadly cons istent with first principle s 
models of giant planet formation (iBromlev & Kenvoiil 12010*). 
which suggest that radial migration and dynamical instability 
ought to be common, even prior to the dispersal of the gas 
disk. The specific separation of planets that we have assumed 
matches that expected if multiple giant planets interact with the 
gas disk such that they end up trapped in mean-motion res- 
onances (Snellgrove et al. 2001; Lee & Peale 2002; Kiev et aO 
2004), provided that the majority of these resonances are broken 
before or shortly after the dispersal of the gas disk. Under some 
conditions fluctuating gravitational perturba tions from disk tur- 
bulence may be able to break resonances (lAdams et al l l2008t 
Lecoanet et al. 2009), though the ubiquity of this process re- 
mains unclear since the true strength of turbulence within disks 
cannot yet be reliably determined. 
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In paper 2 we test the effect of the mass distribution of giant 
planets, including systems with much lower-mass giant planets 
that can undergo Nice model-like instabilities and systems with 
equal-mass giant planet s that produce the stro ngest instabilities 
for a given planet mass ([Raymond et alJlToid) . 

2.1.3. Outer Planetesimal Belt 

In all cases, the outer planetesimal disk contains 50 M© spread 
over a 10 AU-wide annulus. The inner edge of the annulus is 
chosen to be 4 linear Hill radii exterior to the outermost giant 
planet, such that the innermost planetesimals are "Hill stable" at 
the start of the simulation. This disk's mass and mass distribution 
is comparable to that invoked by recent models of the evolution 
of the Solar System's giant pl anets (th e "Nice model") to explain 
the la te heavy bombardment (iTsiganis et al. 2005; Gomes et alj 
l2005h . Note that, for numerical reasons and also to account for 
the much larger total mass, the mass of a "cometary" planetesi- 
mal particle in the outer system is roughly an order of magnitude 
larger than an "asteroidal" planetesimal in the inner system, al- 
though in both cases each planetesimal particle is assumed to 
represent an ensemble of smaller bodies (see below). 

In paper 2 we test the effect of varying certain properties of 
the outer planetesimal disk mass, including the width of the outer 
planetesimal disk and the presence of ~ Earth-mass "seeds" in 
the disk. 



2.2. Integration Method 

Our simulations were run using the hybrid in t egrato r 
in the Mercury integration package (Cham bersI Il999h . 
which combines the speed of a symplectic mapping 
scheme dWisdom & HolmanI Il99lh while particles are well- 
separated with the Bulirsch-Stoer method during closer 
encounters. We used a timestep of 6 days for each of our sim- 
ulations. Particles were removed from a simulation when they 
attained perihelion distances of less than 0.2 AU, at which point 
they were assumed to collide with the star, or if they reached 
aphelion at more than 100 AU, when they were assumed to 
be ejected. Collisions were treated an inelastic mergers. We 
performed extensive numerical tests to validate our choice of 
timestep as well as the threshold integration error in orbital 
energy above which simulations were rejected as unreliable. 
These tests are described in Appendix A. 



2.3. Debris Disk Modeling 

To calculate the dust flux in a planetary system fro m a simula- 
tion w e follow the procedure outlined in section 2 of lBooth et al.l 
( |2009|) with only a few small modifications. This approach 
makes the assumption that planetesimal particles act as tracers 
of a coUisional populations of small bodies and adopts a sim- 
ple model for the evolution of the population. The properties of 
this collisional population are drawn from previous studies that 
fit models of the collisional evolution of planetesimal belts to 
the statistics for the evolution of debris disks (Dominik & D eciiJ 
,2003; Krivov et al. 2005, 2006; Wvat t et al.ll2007bl: iLohne et^ 
2008. .Wvatt .2008: .Krivov ,2010: Ka'ins et alJ l201 ih . Although 
the parameters in this simple model have a precise physical 
meaning for the model as it is presented below, in practice they 
represent "effective" parameters that determine the mass evolu- 
tion and are calibrated to match observations. For example, al- 
though we adopt the classic single power law size distribution 



for a collisional cascade (lDohnanvilll969HWilliams & Wetherilll 
11994 . the true size distribution of a collisional population 
is certainly more complicated (O'Brien & Greenberg 20051 
iKobavashi & Tanakall2010h . We also assume a constant value 
for the impact energy required to catastrophically disrupt 
an object of size D that i s not realistic because g* is undoubt- 
edly a function of D (e.g.. lBenz & AsDhaug'1999). However, this 
type of simple model does a reasonably good job at matching 
more detailed calculations of debris disk evolution that include 
a size dependent an d (consequently) a multiphase size dis- 
tribut i on (see Fig. 1 1 of Lohne et ai]|2008l; iKenyon & BromlevI 
20081 120101; iKrivovl 12010'). Furthermore the parameters of the 
mod el have been tuned to match the observ ed debris disk evolu- 
tion dWvatt et al.ll2007btlKains etal]|201 Ih . 

We divide our planetesimal populations into two compo- 
nents: asteroidal and cometary: asteroids are simply planetesi- 
mals interior to the giant planets' initial orbits and comets are 
exterior We assume that each planetesimal population is in a 
collision-dominated regime from the start of the simulation such 
that its differential size distribution is represented by n{D) oc 
£)2-3^rf^ where D is the ob ject diameter an d - 11/6 for an in- 
finite collisional cascade jPohnanvill 9691 ; i Williams & Wetherilll 
1994) . This distribution spans from the smallest assumed par- 
ticles Dhi = 2.2 fim up to the largest Dc = 2000 km. Smaller 
particles than Dhi are assumed to be blown out of the system by 
radiation pressure on short eno ugh tirnescal es to not contribute 
to the size distribution (see, e.g.. IWv atll2005h . The largest bodies 
in the distribution were chosen to match the largest known ob- 
jects in the Kuiper belt, assuming that the primitive Kuiper belt 
represents a proxy for the accretion that would have occurred in 
planetesimal disks of this mass. If Dc were significantly smaller 
then the collisional timescales would be shorter and the colli- 
sional evolution would proceed more quickly. 

The surface area of a population of bodies with a collisional 
size distribution can be easily calculated. Using our chosen val- 
ues for Dhi, Dc and q^i, and assuming a physical density for all 
objects of 1 gcm"^, the total cross-sectional area in a given ra- 
dial bin is related to the total mass in that bin by: 



M(R) 



= 0.19AU^M„', 



(3) 



where o-(R) is in AU^ and M(R) is in M®. 

For each population of bodies (asteroidal and cometary), 
we calculate the collisional timescale tc, which is a function 
of the parti cle size D as we l l as th e orbital properties of the 
population (IWvatt et al.lll999l[2007ah . This represents the mean 
timescale between collisions that are violent enough to destroy 
bodies of a given size at a mean orbital distance Rm with an an- 
nular width dr. 



tc(D) = 



'2[l + 1.25(e//)2] 



-l/2\ 



where 



(Ttot 



cr(R) 
M(R)' 



(4) 



(5) 



Here tc is in years, is the stellar mass in solar masses and 
cr,„, is the total surface area in AU^. The orbital characteristics 
of the planetesimal population are represented by their mean ec- 
centricities e and inclinations / in radians. The factor fcc(D) rep- 
resents the fraction of the total size distribution that could cause 
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a catastrophic colUsion with a particle of size D, and for our as- 
sumptions can be expressed as 



fcciD) = 



3qd-5 



D 



3qd-5 



+ 2D 



Dcc(D) 



d: 



+ D 



2Dcc(D) 



3qd - 3 



(6) 



where DcAD) - XcD where XcD > Du and Dcc(D) - Dti other- 
wise, and 



Xr = 1.3 X 10" 



-1 n1/3 



1.25e2 + /2 j 



(7) 



is the impact energy required to catastrophically disrupt and 
disperse a body of size D. Here we assume that Q*^ - 200 J kg ' 
and that it does not vary with D. As discussed above, this as- 
sumption is not realistic in terms of the collisional dynam- 
ics (.Benz & Asphaug..l999.) . but for our purposes 2* represents 
an effective strength that determines the dust production and 
mass loss rate from the planetesimal belt, and using a constant 
value allows for a goo d fit t o the statistics of debris disks around 
A stars ( Wva tt et al.ii2007bl) . In addition, it may not be reason- 
able to assume the same 2* for asteroidal and cometary plan- 
etesimals given that their compositions are likely to be different. 
However, the typical collisional timescales for the largest (2000 
km) bodies are very short, just a fewxlO'* years, for the aster- 
oidal planetesimal population (compared with 1-3x10** years 
for the cometary planetesimal population). Thus, asteroidal dust 
is severely depleted within the first 0.1-1 Myr of each system's 
evolution and our choice of 2* has virtually no effect on the 
results. 

For a given simulation timestep, we calculate the spectral 
energy distribution (SED) of the system dust as follows (again, 
as in Booth et al. 2009 with only small modifications): 

- Calculate tc for the largest bodies in both asteroidal and 
cometary planetesimal populations. We then artificially de- 
crease the mass of the planetesimals in each population by 
a factor of [1 + t/tc(Dc)]^^ to account for collisional mass 
loss. This represents a very slow decrease for the comets but 
a significant mass loss among the asteroidsQ. 

- Divide the radial domain into A^^,,,, radial bins spaced log- 
arithmically between 0.2 and 100 AU, which are the inner 
and outer boundaries of our simulation. We tested a range of 



^ Our procedure of decreasing the particle mass according to the col- 
lisional timescale is not completely self-consistent as our planetesimal 
mass did not change during the simulation itself. Given the long colli- 
sional timescales in the outer disk, this assumption has little to no effect 
on the outer planetary systems. This assumption also has little effect on 
cases in which the giant planets were unstable quickly. However, for 
stable simulations or simulations with delayed instabilities, this means 
that, in the inner disk, our dynamics is not completely true to our cal- 
culated dust flux. In other words, the amount of damping provided to 
growing embryos by planetesimals (via viscous stirring and dynamical 
friction) should realistically have been lower if we had accounted for 
collisional evolution of planetesimals. However, given that we compare 
our simulations with observed debris disks that are generally far older 
than the 10-100 Myr timescale for terrestrial planet formation, this does 
not affect our results. For a car eful treatment of dust produ ction during 
terrestrial planet formation, see lKenvon & BromlevI (|2004^ 



Nbin values and found good convergence with Nbi„ > 30 - 50 
and so we use A^;,,-,, - 100 to be conservative. 
Calculate the total planetesimal mass in each radial bin for 
asteroids and comets by sampling the orbit of each plan- 
etesimal Nsamp times along its orbit at intervals that are 
equally spaced in time (i.e., in mean anomaly). We calcu- 
late A^, 



samp 



1 +e/eii„i,, where e//,„/, represents the minimum 
statistical eccentricity needed to cross between radial bins: 
eiimii ~ 0.06 in our case but we use half of that value. Thus, 
we sample circular orbits sparsely but more eccentric orbits 
up to 30 times per orbit to allow the dust to be spread over a 
range of orbital distances. 

For each radial bin and for asteroids and comets, calculate 
the blackbody temperature 



Thb = 278.3 



1/4 



llAU/ 



-1/2 



(8) 



where is the stellar luminosity, R is the radial distance, 
and Tbh is in Kelvin. Then, assuming that all objects radiate 
as blackbodies, we calculate the flux density (in Janskys) Fy 
seen by an observer at distance d: 



Fv = ^2.35xlO""cr(7;): 



BJA,Ti,b(R))\—j X-/, 



-2 



(9) 



where By is the Planck function in units of Jy sr"' and X^ 
is a factor derived by Wyatt et al. (2007b) to account for 
a decrease in emission beyond ~ 200/im needed to match 
observed sub-mm measurements of debris disks: X^ - 1 for 
A < 210;um and X^ = A/21Q for longer wavelengths. 
- Sum the contributions from each radial bin over the whole 
SED, and then sum the asteroidal and cometary contribu- 
tions. We consider the wavelength range from 1 to 1000 yum. 

Using this method, we calculate the SED of each simulation 
from the orbits of all planetesimal particles at each simulation 
timestep. To make the SED useful for comparison with observa- 
tions, the most convenient quantity is the dust flux relative to the 
stellar flux. The stellar flux Fyi, in Jy can be calculated as 



/^. = i.77z..(i,r.)(^)r-(A 



-2 



(10) 



where By is again the Planck function in units of Jy sr ' and 
is the stellar effective temperature, which we take to be 5777 K 
for all cases. Armed with the stellar flux as well as the dust flux, 
we can isolate the dust-to-stellar flux ratio at any wavelength of 
interest such as those observed with Spitzer or Herschel. 

Our simulations ran for 100-200 Myr, but we want to com- 
pare fluxes with stars at a range of ages. Thus, we need to ex- 
trapolate the expected dust fluxes into the future. We do this by 



^ Small grains with sizes of 2.2/jm do not actually radiate as black- 
bodies. Rather, the effective temperature of small grains is likely to be 
higher than the blackbody temperature. The co nsequence of this in the 
context of our m odel is discussed extensively in lBonsor & Wvat3 d2010h 
and Kains et al. i201 ll ). To summarize, the actual flux may be slightly 
higher or lower than that predicted with our assumptions, since on the 
one hand the higher temperature means that we underestimate the dust 
flux for a given dust mass, however this is compensated to some extent 
by the fact that we also overestimate the evolutionary timescale (since 
the evolutionary timescale for a disk of observed temperature is set ob- 
servationally). 
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making the assumption that there is no more significant dynami- 
cal evolution of the system, i.e. that the orbital characteristics of 
the planetesimals are not going to change in the future. Clearly, 
this assumption does not account for pun ctual events like the late 
heavy bombardment (iGomes et al ] l200l or very late dynamical 
instabilities. For a given time t after the end of the simulation, 
we decrease the planetesimal mass by a simple factor related 
to the coUisional timescale at the end of the simulation, i.e., as 
(1 + f/fc) In terms of the global analysis of our sample of sim- 
ulations, we retain snapshots of each simulation at wavelengths 
of 1, 3, 5, 15, 20, 25, 50, 70, 100, 160, 250, 350, 500 and 850 ^lm 
at time intervals of 1, 3, 10, 30, 100, 300, 1000, and 3000 Myr 

Our dust flux calculations compare reasonably well with 
observations and other models. We perform this comparison 
using systems that are dynamically calm (referred to as "sta- 
ble systems" later in the paper), in which the giant planets re- 
main on stable orbits throughout the simulation such that the 
outer planetesimal disk remains largely intact. In these stable 
systems, our dust flux calculations yield typical values for the 
dust-to-stellar-flux ratio FjFstar of 0.1-0.5 at 25//m and 10- 
35 at 70yum after 1 Gyr of simulated dynamical and calcu- 
lated coUisional evolution. These values are broadly consistent 
with the fluxes detected around solar-type star s with observed 
exces s es at 24 and 70^ m {Habing et al. 2001; Beichman et al 
12006; 'Moor et al." 2006; Trinina et al. 2008; Hiflenbrand et al 
[2008; Carpenter et al. 2009; Caspar etal. 2009), although our 
stable simulations yield very few systems with F/FstarPOjjm) x 
1-10, probably because of the relatively large masses in 
our outer planetesimal disks (note that our unstable systems 
can yield those flux levels). Compared with the more sophis- 
tic ated_jnodels_of_du^t prod^ planetary accretion 
of lKenvon&BromlevI ( I2008L 1201 Ol) . our calculated fluxes are 
larger by a factor of a few, notably at 70yum. Our fluxes are only 
a factor of 2-3 larger than those calculated bv iKennedv & Wvatti 
(I2OIOI) . 
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Fig. 3. The eccentricity (top panel) and inclination (bottom 
panel) evolution of the two surviving giant planets in the ex- 
ample simulation. The three curves in the right panel show the 
inclination of each planet relative to the starting plane of the sys- 
tem - presumably the stellar equatorial plane - in black and light 
grey and their mutual inclination, which is larger than each of 
their individual inclinations, in dark grey. 



3. An example simulation 

In this section we explore the dynamics and dust properties of 
a single simulation. In following sections we will consider the 
outcome of an ensemble of many simulations, so this is an op- 
portunity to inspect the detailed evolution of a particularly inter- 
esting system. The initial conditions for the chosen simulation 
are shown in Fig. [T] the giant planet masses were 1.46 Mj (at 
5.2 AU), 0.75 Mj (7.9 AU) and 0.55 Mj (11.3 AU). 

Figure[2]shows the dynamical evolution of the system, which 
became unstable after 21 million years. Before the instabil- 
ity, the system evolved in the expected quiescent fashion. In 
the inner disk, embryos accreted planetesimals and other em- 
bryos from the inside-out. Embryos' eccentricities and inclina- 
tions were kept small by the planetesimals via dynamical fric- 
tion and viscous stirring. After 21 million years there were sev- 
eral almost fully-grown terrestrial planets, including three plan- 
ets more massive than 0.6 at 0.61, 0.93 and 1.34 AU and a 
dozen other embryos of 0.06-0.3 M® extending out to 2-3 AU. 
During this phase of quiescent accretion, the giant planets' or- 
bits remained almost perfectly circular, with eccentricities of less 
than 0.05 for each planet (less than 0.03 for the innermost gas 
giant). Planetesimal scattering caused a slow drift in the giant 
planets' semimajor axes, of 0.06 AU for the innermost (and most 
massive) giant planet, and less than 0.01 AU for the other two 
giants. During this time, the outer planetesimal disk was slowly 
sculpted by the giant planets. The inner edge of the planetesimal 
disk was eroded, and strong resonances with the outermost giant 



planet (low-order mean motion resonances) were slowly cleared 
out. At the time of the instability, roughly three quarters of the 
initial outer planetesimal disk was intact. 

The instability began with a rapid eccentricity increase in 
the eccentricity of the two outer giant planets and was followed 
by a close encounter between those two. That set off a series 
of close encounters between all three giant planets that lasted 
84,000 years, involved 35 planet-planet scattering events, and 
culminated with the ejection of the outermost giant planet. The 
behavior of the simulation was characteristic in terms of dy- 
namical instabilities in that the least massive giant planet was 
ejected and the surviving planets were segregated by mass, with 
the most massive one closest to the star (a t 3.65 AU) and the less 
massive one fa rther out (at 36.6 AU, e.g. lChatteriee et al1[26()8l: 
[Raymond et aL.201Q) . 

During and immediately after the giant planet instability, the 
inner disk was almost entirely thrown into the star This hap- 
pened by a rapid increase in embryos' and planetesimals' ec- 
centricities mainly by secular forcing from the scattered giant 
planetsQ Among the 3.9 of terrestrial material that was de- 
stroyed were 13 embryos including a 0.95 Me at 0.6 AU, a 0.62 

In the simulation, any body that came within 0.2 AU of the star is 
considered to have collided with it. It is conceivable that in some cases, 
star-grazing embryos could have their orbits shrunk and re- circularized 
by tidal effects, altho ugh the efficiency of th is process is uncertain and 
probably quite small jRavmond et alj2008bh . 
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Fig. 2. Evolution of a system with a relatively late (21 Myr) instability among the giant planets. Each panel shows a snap- 
shot in time of orbital eccentricity vs. semi-major axis for all particles; vertical bars denote sin(/) for terrestrial bodies with 
Mp > 0.2 Me and / > 10°. The particle size is proportional to the mass'^^, but giant plane ts (large black circles) are 
not on this scale. Colors denote water content, assuming a Solar System-like initial distribution (iRavmond et al.l 12004). The 
surviving terrestrial planet has a mass of 0.72 M®, a stable orbit within the habitable zone (semimajor axis of 0.96 AU), 
and a high eccentricity and inclination (and large oscillations in these quantities). A movie of this simulation is available at 
[http://www.obs.u-b ordeauxl.fr/e3arths/raymond/scatterSED.mpg 



Me planet in the habitable zone at 1.34 AU, and five other em- 
bryos larger than 0.15 Me- All but one of the embryos - the 
sole survivor - collided with the star within a few hundred thou- 
sand years of the start of the instability. The rocky planetesimals 
were all destroyed within 1 Myr The outer planetesimal disk 
was completely destabilized by the instability, as the two lower- 
mass giant planets scattered each other to tens of AU. All plan- 
etesimals were either ejected from the system (about 85%) or 
collided with the star (15%). The vast majority were destroyed 
within 1 Myr of the instability - all but three within 5 Myr - but 
the last planetesimal took almost 25 million years to be ejected. 



Three planets - two gas giants and one terrestrial planet - 
survived to the end of the simulation, all on excited, eccentric 
and inclined orbits. The long-term dynamics of the two surviv- 
ing giant planets is shown in Figure [3] The eccentricities of the 
surviving giant planets are considerable as are their orbital in- 
clinations with respect to the starting plane. The two giant plan- 
ets have a mutual orbital inclination of ~ 36°. Although this is 
less than the formal limit of 39.2° required for the Kozai effect 
to take place (Kozai 1962; Takeda & Rasio 2005), out of phase, 
Kozai-like oscillations are evident in the planets' eccentricities 
and inclinations. 
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Fig. 4. The eccentricity (top panel) and inclination (bottom 
panel) of the surviving terrestrial planet in the example simu- 
lation (mass of 0.72 Me, semimajor axis of 0.93 AU). The two 
curves in the right panel show the inclination relative to the start- 
ing plane of the system in grey and the mutual inclination with 
respect to the innermost giant planet in black. 



The terrestrial planet's final mass is 0.72 M© and its semi- 
major axis of 0.935 AU places it in the cir cumstellar habitable 
zone (iKasting et alJll993tlSelsis et alJl2007 l). However, Figure g] 
shows that the planet's orbit is strongly perturbed. Its eccentiic- 
ity oscillates between 0.4 and 0.7 and its inclination between 
almost zero and more than 60° on a ~ 10,000 year timescale 
(although there other longer secular frequencies in the oscilla- 
tions). On Myr timescale the orbit has a minimum and maxi- 
mum eccentiicity of 0.39 and 0.73 and a minimum and maxi- 
mum inclination of 0.03° and 63.3°. Given its large eccentiicity, 
one would e xpect this planet's climate t o be highly variable d ur- 
ing the year dWilliams & Pollardll2002HDressing et aDl2010h . In 
addition, the changes in both its eccentricity and inclination - 
equivalent to changes in obliquity for a fix ed sp in axis - woul d 
cause variations on the secular timescales (ISpiegel etliDl2010h . 
Given that the orbit-averaged stellar flux increases with the or- 

bital eccentricity (as (1 - e l ) and that the planet's closest 
approach to the star is only 0.25 AU, this planet may not be 
habitable during its high-eccentricity phases. However, more de- 
tailed modeling of such planets is beyond the scope of the pa- 
per and the reader is referred to recent climate modeling papers 
that include the effects of var ying the eccentricity and obliq- 
uitv Jwil liams & Poflard 2001 120031: ISoiegel et al.ll2009l l20Ta 
[Dressing et al.ii201 0). 

The evolution of the dust brightness follows from the dynam- 
ical and collisional evolution of the planetesimals in the system. 
Figure |5] shows the resulting spectral energy distribution at dif- 
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Fig. 5. Evolution of the spectral energy distribution (SED) of 
the simulation from Fig. |2] Each curve shows the SED during 
a given simulation snapshot. The instability occurred at 2 1 Myr, 
and the SED evolved dramatically in the immediate aftermath, 
as icy planetesimals were scattered onto high-eccentricity orbits 
- thereby producing transient hot dust - and then ejected from 
the system. The dashed line represents the stellar photosphere. 
The system is viewed at a distance of 10 pc. 



ferent snapshots in time. The collisional timescale for the largest 
planetesimals (2000 km) in the terrestrial planet forming region 
is only tcoii ~ 10'* years, so the planetesimals in the inner disk 
are ground into hot dust within just a few million years. Thus, the 
total dust brightness has dropped sharply at all wavelengths by 
the 20 million year snapshot, most dramatically at wavelengths 
shorter than roughly 50 fim. For the rest of the simulation, the 
primary source of dust is the outer planetesimal belt for which 
tcoii ^10^ years because the inner disk planetesimals have been 
ground awayO 

When the instability occurs, the spectral energy distribution 
changes dramatically and quickly (Fig. |5]l. Comets from the 
outer planetesimal disk are scattered onto highly eccentric or- 
bits. Those that enter the inner Solar System can cause bombard- 
ments on the terrestrial planets akin to or often far more intense 
than the Solar System's late heavy bombardment (iGomes et al.l 
120051: IStrom et al.ll2005? ). although given the much earlier tim- 
ing of most instabiUties, these bombardments could act to seed 
the terrestrial planet-forming region rather than impact fully- 
formed planets. The decrease in the comets' perihelia introduces 
a large amount of warm dust into the system which lasts for 
the duration of the bombardment and leads to a spike in bright- 
ness at near- to mid-infrared wavelengths, shown in Figure )6]for 
A ^ 5,10,25,70,16 and 500yum (see also iBoofli etal J 120091: 
iNesvornv et al.ll20Toh . At shorter wavelengths the flux is more 
jagged than at longer wavelengths because the short wavelength 



' We do not consider regeneration of small bodies via giant impacts 
but we note that the debris from giant collisions could cause short-lived 
spikes in the dust brightness, in particular at mid- to near-infrared wave- 
lengths jS tem 1994; Grogan et al. 2001: iKenvon & BromlevI |2005| ; 
iLisse et al.j,200 9). These peaks in brightness from collisions can only 
occur within a few AU because a very massive collision is needed to 
cause substantial brightening. 
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Fig. 6. The ratio of the dust-to-stellar flux as a function of time for six different wavelengths from 5 to 500 microns, including a 
zoom during the instability - note that each main plot is on a log scale and each zoom-in is on a linear scale. The rough observational 
limits of the MIPS instrument on NASA's Spitzer Space Telescope are shown for 25 and 70 fim with the dashed line iTrilling et alj 
0200 8). All planetesimal particles were destroyed as of 45 Myr via either collision or ejection so there is no more dust in the system 
after that point. 



flux is entirely due to hot dust that is produced sporadically 
by individual planetesimals entering the inner planetary sys- 
tem. In contrast, the long wavelength flux combines the flux 
from a much larger range of dust temperatures and therefore 
includes a much larger number of particles. The spike in flux 
associated with the instability is strong for wavelengths shorter 
than ~ SOfim, but at longer wavelengths the spike is weak or 
absent. As objects are dynamically removed from the system 
the system's brightness drops precipitously and out of the de- 
tectable range in the few million years following the instabil- 
ity. We note that a more realistic model of dust production by 
new comets suggests that the mid-infrared peak during the bom- 
bard ment in our model ma y be underestimated by a factor of a 
few (iNesvornv et al.ll2010h . 



4. Results for the ensemble of simulations 

In this section we study the outcome of our ensemble of simu- 
lations (called mixed in paper 2). This set is of particular inter- 
est because the surviving giant planets match the observed ec- 
centricity distribution without any fine-tuning (Raymond et al. 
l2009al l2010() . We explore the formation process of terres- 
trial planets, the orbital characteristics of terrestrial planets that 
formed, giant planet dynamics, dust production from planetesi- 
mals, and correlations between these. We also compare the sim- 
ulations with the observed extra-solar giant planets and debris 
disks. In paper 2 we explore the effect of a variety of parameters 
on the process. 



4. 1 . Giant Planet Instabilities 

Of the 152 simulations that ran for at least 100 Myr and met 
our energy conservation cutoff (see Appendix A), the giant plan- 
ets were unstable in 96 systems (63%). The instability times 
ranged from 247 years to 180 Myr after the start of the simu- 
lation with a median of 91,600 years. This is encouraging be- 
cause it is close to the value expected based on our initial gi- 
ant planet separations of 4 - 5RH,m dMarzari & WeidenschillingI 
120021: IChatteriee et al.ll2008l) . This means that our assumption of 
no disk gas at the start of the instability is marginally acceptable, 
because the typical timescale for the gas dissipation is thought 
to be ~ 10^ years ( Wolk & Walter 1996; Ercolano et al. 201%. 
However, in many cases instabilities are likely to have occurred 
in the presence of a residual gas disk. To test this assumption, 
we ran an additional set of simulations that include a simple pre- 
scription for the effects of gas damping that are presented in pa- 
per 2 (and show no significant changes from the gas-free simu- 
lations). We note that later instabilities can certainly occur (e.g., 
the late heavy bombardment) but we expect the number of in- 
stabilities to decrease in time, and the required computing time 
made it impractical to search for instabilities beyond 100-200 
Myr 

Although 75% of instabilities occur within 1 Myr, there is 
a tail that extends to longer timescales. One in six instabilities 
(16 out of 96) occurred after 10 Myr, one in ten (10/96) after 30 
Myr, and one in fifty (2/96) after 100 Myr Given that our simu- 
lations only lasted 100-200 Myr, we undersampled the fraction 
of unstable systems in the 100-200 Myr timespan and we expect 
that even later instabilities should certainly occur in a fraction of 
systems due to long-term chaotic diffusion. The timing of the in- 
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Fig. 7. Eccentricity distributions of the surviving giant (solid 
black line) and terrestrial (dashed black line) planets in the un- 
stable systems. The grey curve represents the known extra-solar 
planets beyond 0.2 AU - this helps to exclude planets whose 
orbits have been tidally circularized as well as low-mass planets. 



Fig. 8. The total mass in surviving terrestrial planets as a func- 
tion of the minimum perihelion distance of a giant planet during 
the simulation. Black dots represent systems in which the gi- 
ant planets were dynamically unstable and grey dots are systems 
that were stable. Filled black dots are systems in which a single 
terrestrial planet formed. 



stability is important in terms of the sizes of objects in the inner 
disk; instabilities later than the typical terrestrial planet forma- 
tion time of 10-100 Myr may destroy fully-formed Earth-sized 
planets - sometimes with appreciable water contents - rather 
than embryos and planetesimals. 

As shown in Figure |7] the surviving giants in the unstable 
simulations provide a quantitative match to the observed extra- 
solar giant planets (p = 0.49 from a Kolmogorov-Smirnov test, 
consistent with previous work with very similar initial condi- 
tions but much larger samples; iRavmond et allBOOSal l2009bla[ 
12010 ). Thus, one might imagine that the unstable systems rep- 
resent the appropriate subsample of simulations that we should 
use to represent the known exoplanet systems, and that our stable 
simulations are unrealistic in that they somehow la ck a trigger to 
make them unstable (e.g.. iMalmberg et al.ll2010l) . However, as 
we discuss in section 5.1, the observed exoplanets do appear to 
require a contribution from dynamically stable systems. We note 
that additional constraints exist in the exoplanet sample (e.g., the 
mass-eccentricity correlation). Combinations of different sets of 
simulations can match all of the observed characteristics of the 
exoplanet distribution, and we perform this exercise in paper 
2 (see also section 5 in Rav mond et alJl20Toh . 

A prediction of the planet-planet scattering model is 
that most planetary systems should contain multiple gi- 
ant planets, i.e., additional giant planet s should exist exte- 
rior to most of the kno\yn one s (e.g., iRasio & FordI Il996t 
iMarzari & Weidenschillin jl2002l) . Likewise, the vast majority 
of unstable simulations in our sample (85/96 = 89%) contain 
multiple giant planets; the remaining 11% contain just a single 
surviving giant planet. The outer giant planet is typically 5-10 
AU more distant than the inner giant planet (with a tail to > 30 
AU), so long-duration observations are needed to follow up the 
known giant exoplanets. The distribution of separations between 
planets in scattered systems may actually provide constraints on 
their initial mass distribution because the surviving planets tend 
to be more widely-sepa rated if they are equal-m ass than if their 
masses differ markedly (iRavmond et al.ll2b09bl) . 



4.2. Terrestrial Planet Formation 

The number and spacing of terrestrial planets that form is gov- 
erned by th e eccentriciti es of the planetary embryos from which 
they form (iLevison & A gnor 2003). These eccentricities are a 
result of gravitational forcing both from nearby embryos and 
the giant planets. There are two differences in terrestrial planet 
formation between stable and unstable systems: 1) perturba- 
tions during the instability can play a major role in shaping 
the embryo distribution; and 2) the dynamical state of the gi- 
ant planets after an instability is generally more excited than be- 
fore the instability. In some unlucky cases the surviving giant 
planets provide an accretion-friendly environment that is empty 
because the instability has already removed all rocky bodies. 
Past simulations of terrestrial planet formation with giant planets 
on excited orbits have all neglected the planet-planet scattering 
phase during which the giant planets a ctually acquired their ec- 
centricities jCh ambers & Cassen 200^ iLevison & Agnorll2003l 
[Raymond et alJ l2004t (Raymond 2006b . which clearly plays a 
very important role in the dynamics. 

Giant planet perturbations span a continuous range but the 
outcome is quantized into a discrete number of terrestrial planets 
during the accretion process. If perturbations are weak - if the 
giant planets collide rather than scattering (or are dynamically 
stable or low-mass) - then embryos' eccentricities remain small 
and feeding zones narrow and several terrestrial planets form. 
For stronger giant planet perturbations, feeding zones widen and 
fewer terrestrial planets form, although the total mass in planets 
tends to decrease because stronger perturbations imply that the 
giant planets were scattered closer to the terrestrial planet region 
so more embryos end up on unstable orbits. In systems where 
embryos' radial excursions are comparable to the radial extent of 
the surviving disk only one planet forms, usually on an excited 
orbit. In the simulation from Fig. |2] the lone terrestrial planet 
did not accrete from a disk of excited embryos but rather was 
the only planet to survive the instability. Perturbations during, 
not after, the instability determined the outcome. 

The strength of the giant planet perturbations is directly re- 
lated to the smallest giant planet perihelion distance qcp,mm- 
Figure |8] shows that the efficiency of terrestrial planet formation 
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Fig. 9. Distribution of the number of surviving terrestrial planets 
(defined as Mp > 0.05 M©) in our fiducial mixed set of simula- 
tions. The grey shaded histogram shows the unstable simulations 
and the dashed line shows the stable simulations. 



Fig. 10. The peak-to-peak oscillation amplitudes in eccentricity 
e and inclination / for the surviving terrestrial planets in simula- 
tions with stable (grey) and unstable (black) giant planets. The 
single-terrestrial planet systems are shown with filled black cir- 
cles, and Earth is the grey star. 



is directly related to qcp,min- This is not surprising: planets with 
smaller qcp.min have higher eccentricities (since all giant planets 
start with the same range of semimajor axes) and have therefore 
undergone more violent scattering events. Every system with a 
giant planet scattered to qcp,mm < 1 .3 AU destroyed all terres- 
trial material in the system, as did so me simulations out t o 3 AU . 
This is consistent with the results of I Veras & Armitage' (l2006h . 
who found a similar scaling between qGP,mm and the survival of 
test particles in the inner disk. As expected, systems that form a 
single terrestrial planet are those with giant planet perturbations 
that are almost strong enough to completely empty the terrestrial 
material from the system. In these cases the embryos' eccentric- 
ities are large and, as we will see below, the single planet that 
survives maintains an excited orbit. 

4.3. Characteristics of surviving terrestrial planets 

All terrestrial material was destroyed in 41 out of 96 unsta- 
ble simulations (43%; Figure|9]l. Likewise, 22 unstable systems 
(23%) formed a single terrestrial planet, 16 formed two terres- 
trial planets (17%), and and the remaining 17 systems (18%) 
formed three or more terrestrial planets. Among the 56 stable 
simulations, only 7 (12.5%) formed two planets, and the remain- 
ing 87.5% of simulations formed three or more planets. 

As a population, the surviving terrestrial planets have smaller 
eccentricities than the giant planets (Figure |7}; the median ec- 
centricity is 0.10 for terrestrials and 0.21 for giants. The terres- 
trial planets that formed in systems with unstable giant planets 
have only slightly more eccentric and inclined orbits than the ter- 
restrial planets that formed in stable systems. The median e and 
; were 0.10 and 5.1° for the unstable systems and 0.08 and 2.5° 
for the stable systems, respectively. The single-terrestrial planet 
systems were the most excited, with median e of 0.14 and / of 
8.7°. If we neglect the single-terrestrial planet systems, the ec- 
centricities and inclinations of the terrestrialplanets that formed 
in stable and unstable systems is very closeQ 



* We expect that numerical resolution should play a role here. Given 
that they contain only 500 planetesimal particles, our si mulations can- 
not fully resolve dynamical friction at late times (e.g. lO'Brien et al.l 
l2006l : iRavmond et al.ll2b06bh . We expect that eccentricities and incli- 



Although their time-averaged orbits are similar, terrestrial 
planets in unstable systems undergo much stronger orbital os- 
cillations (Figure [Toll. The median peak-to-peak eccentricity and 
inclination oscillation amplitudes are 0.11 and 5° for the unsta- 
ble systems and 0.043 and 2.8° for the stable systems, respec- 
tively. Again, the single-terrestrial planets are the most excited 
and undergo the largest oscillations, with median e and / ampli- 
tudes of 0.21 and 11.9°. The climates of the single planets are 
likely to vary dramaticall y on the secular timescale of 10^ - 10* 
years (ISpiegel et al.ll2010l) . 

4.4. Correlations with observable quantities: giant planets 
and debris disks 

Observations of massive exoplanets can only very roughly di- 
agnose the outcome of terrestrial accretion. Figure [TT] shows 
a strong negative correlation between the efficiency of terres- 
trial planet formation and the eccentricity of the innermost giant 
planet eg. For stable systems gg tends to be very small (median 
eg = 0.008) while for the unstable systems eg spans a very wide 
range, from <0.01 to 0.8 (median eg - 0.21). The correlation 
between the total terrestrial planet mass and eg is expected but 
the range in terrestrial planet mass for different systems with a 
similar eg shows the importance of the orbital history. 

The orbits of surviving giant planets retain a memory of the 
strength of the instability, or lack thereof. Figure [T2]breaks down 
the giant planet eccentricity distribution by outcome in terms of 
the number of terrestrial planets that formed. Multiple terrestrial 
planets form preferentially in systems with low giant planet ec- 
centricities, because these represent very weak instabilities. By 
the same argument, highly eccentric giant planets tend to destroy 
all terrestrial material in their systems. The intermediate regime 
is represented by the single-terrestrial planet systems; for these 
systems eg is typically in the range 0.1-0.3 (median of 0.17). 
These ecc entricities are close to the median of the e xoplanet dis- 
tribution (But ler et all 120061: lUdrv & San"tosll2007l) and there is 
considerable overlap from systems with zero or many terrestrial 
planets. It is therefore very difficult to diagnose a single-planet 

nations are somewhat overestimated in the stable systems, but since the 
instabilities remove most of the planetesimals, not in unstable systems. 
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Fig. 11. The total mass in surviving terrestrial planets as a func- 
tion of the eccentricity of the innermost surviving giant planet. 
The Solar System is shown as the grey star. 




Giant planet eccentricity 

Fig. 12. Distribution of eccentricities of the innermost giant 
planet in simulations which formed zero (blue line), one (green), 
and two or more (red) terrestrial planets. The sum of the three 
distributions (black dashed line) provides a quantitative match 
(p=0.49 from a Kolmogorov-Smirnov test) to the observed exo- 
planets (thick grey line). 



system from a measurement of just the eccentricity of a giant 
exoplanet. 

The outer disk's evolution is also governed by giant planet 
perturbations. Icy planetesimals - whose collisional erosion cre- 
ates debris disks - survive in dynamically calm environments 
where the giant planets were either stable, low-mass, or un- 
derwent a relatively weak instability. Indeed, Figure [13] shows 
a strong anti-correlation between the lQfj.m dust flux and the 
giant planet eccentricity. This arises simply because eccentric 
giant planets have increased apocenter distances and impinge 
on the planetesimal disk, thereby dynamically removing plan- 
etesimals via ejection. In addition, if planetesimals survive on 
highl y-eccentric orbits their collisional lifetimes may be short- 
ened (IWvatt et aLllioTob . 

Given that the efficiency of terrestrial planet formation anti- 
correlates with the giant planet eccentricity (Fig. [TTt and the 
dust flux also anti-correlates with the giant planet eccentric- 
ity. Figure [14] shows that the efficiency of terrestrial planet for- 



Fig. 13. The dust-to-stellar flux ratio at 70 //m after 1 Gyr of dy- 
namical and collisional evolution, plotted as a function of the 
eccentricity of the innermost surviving giant planet. Unstable 
systems are plotted in black and stable systems in grey (sys- 
tems with a single surviving planet are shown with solid circles). 
The pileups very close to the x axis represent systems with vir- 
tually zero 70 //m flux. The dashed line represents an approxi- 
mate threshold a bove which excess e mission was detectable us- 
ing Spitzer data dTrilling et al.ll2008l) . The star shows the esti- 
mated flux from the Kuiper Belt at 1 Gyr, as calculated previ- 
ously ((Booth et al. 2009) based on dyn amical simulations of the 
outer Solar System dGomes et al.ll2005h . 



mation correlates with the dust flux. Stars with bright dust al- 
most all contain terrestrial planets: the median system for which 
F IF star > 10 contains 3.6 in terrestrial planets, and every 
single system contains at least 1.85 Me in terrestrial planets. 
Systems that are extremely bright at long wavelengths should 
therefore be considered prime targets in the search for terrestrial 
planets. 

The debris disk-terrestrial planet correlation seen in Fig. [14] 
exists because the inner and outer planetary system are subject 
to the same dynamical environment: the violent instabilities that 
abort terrestrial planet formation also tend to remove or erode 
their outer planetesimal disks. As was the case for the giant 
planet eccentricities, single-terrestrial planet systems populate 
the intermediate area of the correlation and overlap with sys- 
tems with no planets as well as those with several. The terrestrial 
planet - debris disk correlation is not perfect as there exist "false 
positives" with bright dust emission and no terrestrial planets, 
corresponding to systems with asymmetric, inward-directed gi- 
ant planet instabilities. Conversely, "false negatives" with ter- 
restrial planets but little to no dust are systems that underwent 
asymmetric but outward-directed instabilities. We discuss these 
caveats in detail in paper 2 when we compare our simulations 
with the observed debris disk statistics and use them to derive a 
crude estimate of rjEaith ■ 

The debris disk-terrestrial planet correlation is a function of 
wavelength. Figure [15] shows the correlation at 1 Gyr for six 
wavelengths between 5 and 500 jim. At 5;um there is no cor- 
relation but mainly a scatter plot. The only hint of a correlation 
is a lower envelope of dust fluxes larger than F IF star S 10"^° for 
systems with more than ~ 2 in terrestrial planets. This non- 
correlation comes from the fact that the short-wavelength flux 
can be dominated by either a small amount of hot dust or a large 
amount of colder dust. At 15/im the debris disk-terrestrial planet 
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Fig. 15. The dust-to-stellar flux ratio F/Fstar after 1 Gyr of dynamical and collisional evolution as a function of the total mass in 
terrestrial planet s, for six wavelength s between 5 and 500 microns. The Spitzer observational limits are shown for 25 and 70 jjm with 
the dashed line jTrilUng et al.ll2008l) . In each panel, grey circles represent stable simulations and black circles represent unstable 
simulations. 
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Fig. 14. The dust-to-stellar flux ratio at 70;um after 1 Gyr of dy- 
namical and collisional evolution, plotted as a function of the 
total mass in terrestrial planets. Unstable systems are plotted 
in black and stable systems in grey. The pileups close to the 
X and y axes represent systems with either no terrestrial plan- 
ets or virtually zero 70;um flux. The dashed line represents an 
approximate threshold abo ve which excess e mission was de- 
tectable using Spitzer data (iTrilling et al.ll2008h . The star shows 
the estimated flux from the Kuiper Belt at 1 Gyr, as calculated 
previously (Booth et al. 2009) using model s based on dynami- 
cal simulations of the outer Solar System (iGomes et alJl2005b . 
The Solar System falls into an intermediate region of parame- 
ter space: the giant planets may have been unstable, but only 
weakly so (there is no clear evidence for ejections or planetary 
collisions), while the Kuiper Belt would have been bright for 
several hundred Myr prior to the Late Heavy Bombardment. 



correlation clearly exists, although there are still some outliers 
with large fluxes and small terrestrial planet masses. At longer 
wavelengths these outliers slowly disappear because the signal 
from small amounts of transient hot dust is overwhelmed by the 
large amount of cold dust produced in quiescent outer planetesi- 
mal disks. The debris disk-terrestrial planet correlation is evident 
for all wavelengths longer than ~ 25jum. 

Figure [16] (left panel) shows histograms that represent hori- 
zontal slices through Fig. [15] For each wavelength we chose a 
"detection limit" and tabulated the fraction of systems that were 
deemed detectable as a function of the total terrestrial planet 
mass. These detection limits were chosen for illustrative pur- 
poses and are in many cases unrealistic. For example, the cur- 
rent detection threshold at Spva is probably closer to ~ 10"^, 
but none of our simulations would be detectable at that limit. 
The reader is referred t o other work discussing current dust 
detection lim its at short (lAbsil et al.ll2006t lAke son et al."2009h 
and long (Smi th et al.ll2009l) wavelengths (see also Wyatt 200§. 
The fraction of systems with no terrestrial planets that is de- 
tectable varies significantly between the different wavelengths 
and simply reflects the location of the detection limit with re- 
spect to the general trends in Fig. [15] The detection limits at 
25A tm and 70iuni, whic h correspond to the actual Spitzer lim- 
its (ITrilling et al.ll2008l) . are close to the extremes: at 25;um not a 
single terrestrial planet-free system is detectable while at 70/im 
more than 40% of all systems are detectable. 

At each wavelength, the fraction of systems that is detectable 
increases with the terrestrial planet mass (Fig. [16). At all wave- 
lengths the trend is relatively flat to a certain point where it in- 
creases significantly, by several standard deviations between ad- 
jacent bins. At all wavelengths this transition occurs between 
either the 1-2 M® and 2-3 M® bins or the 2-3 M® and the 
> 3 M® bin. Even at 5;um, which showed no obvious debris disk- 
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Fig. 16. Left: The fraction of systems that would be detectable after 1 Gyr as a function of the total mass in surviving terrestrial 
planets for six different wavelengths between 5 and 500/im. Each curve represents a horizontal slice through Fig [15] at a different 
vertical height. The chosen "detection limits" in FjFsiar were: 10"'^ at Syum, 10""' at ISyum, 0.0 54 at 25iUm, 0.55 at I Qimi, 10 at 
160/im, and 100 at 500/zm. The error bars are 1-cr and were calculated using binomial statistics dBurgasser et al.ll2003h . Note that 
the different curves are offset by up to +0.1 Me for clarity. The bins themselves are at: zero, 0-1 M®, 1 - 2Me, 2 - 3Me, and 
> 3 Me. Right: The fraction of systems with 0.5 Me or more in terrestrial planets as a function of F / F s,ar(70fim) (1 Gyr). Systems 
with F IF star < 10"^ are included in the bin at FjFstar ~ 10"^. The Spitzer detection limit is shown as the dashed line. This represents 
a vertical slice through Fig. [14] 



terrestrial planet correlation in Fig. [15] there is marked jump in 
the fraction of detectable objects in the last bin. This jump is 
also seen at 15;um, although these two shortest wavelengths are 
the only ones for which some systems with more than 3 Me in 
planets are not detectable. Although our chosen detection limits 
at certain wavelengths may be overly optimistic or pessimistic 
(e.g., F/ Fsiari5iJ.m) > 10^^ is likely to be difficult to achieve in 
the near term), this shows that the debris disk-terrestrial planet 
correlation is present for all of these wavelengths. 

The right panel of Fig.[T6]shows the fraction of systems with 
significant mass in terrestrial planets as a function of FjFsiar 
at 70;um. This plot represents a vertical slice through Fig. [14] 
At small dust fluxes, a minority of systems contain terrestrial 
planets - these are referred to in section 5 as "false negatives". 
The fraction of systems with terrestrial planets increases dramat- 
ically beyond the detection limit, and every single one of the 58 
systems with FjFsiar > 10 (including 7 unstable systems) con- 
tains at least 1 .9 Me in terrestrial planets, with a median value of 
3.5 Me. Of the 53 unstable systems above the detection thresh- 
old at IQjjm, 35 (66%;^g g||) contain terrestrial planets. Thus, in 
our simulations debris disks appear to represent signposts for 
terrestrial planets with a confidence of at least 66% at 70/im. 

The anti-correlation between debris disks and eccentric giant 
planets in our simulations also holds across many wavelengths. 
Figure[T7](left panel) shows the fraction of stars that is detectable 
as a function of the giant planet eccentricity for a range of differ- 
ent wavelengths, similar to Fig. [16] Again, the trends are stronger 
for wavelengths of 25;um and longer, but are still clear at shorter 
wavelengths (although the detection thresholds at these wave- 
lengths are certainly much smaller than in reality). Fig. [17] (right 
panel) shows that the fraction of systems containing an eccentric 
giant planet is anti-correlated with the dust emission at IQjjm, 
simply because in our simulations the calmest giant planets de- 



stroy the fewest number of outer planetesimals and therefore 
produce the brightest debris disks. 

The debris disk-terrestrial planet correlation is also a func- 
tion of time. Figure[T8]shows FjFsmr at 70;um vs. the final terres- 
trial planet mass for all simulations at eight snapshots between 
1 Myr and 3 Gyr At 1 Myr all systems are above the detection 
threshold, even those that have akeady undergone violent insta- 
bilities, as the timescale for clearing out planetesimals is gener- 
ally closer to a few to 10 Myr from the time of the instability. 
In a given snapshot, systems that have not yet become unstable 
but that eventually will are those for which the flux remains as 
high as the cluster of stable systems but for which the total ter- 
restrial planet mass is low, indicating a strong future depletion of 
rocky material. After an instability occurs the flux drops (though 
differently at different wavelengths; Fig. [6]l but the total terres- 
trial planet mass, measured at the end of the simulation, does 
not change, such that a given system moves vertically between 
snapshots. In time, instabilities remove systems at low terrestrial 
planet mass and high FjFstar', the last two instabilities occurred 
at 168 Myr (in which one 2.4 M® terrestrial planet survived on a 
highly eccentric orbit) and 180 Myr (in which four fully-grown 
terrestrial planets were destroyed including a 1 .3 Me planet in 
the habitable zone). Beyond the end of our simulations one can 
imagine that a fraction of stable systems could actually become 
unstable and quickly lose a large fraction of their flux (and per- 
haps their terrestrial planets as well). 

Stable systems remain clustered at high fluxes at all times, 
but decreasing in time due to collisional grinding|3 The same 

' Note that one stable system ejected all of its planetesimals. In that 
system the interactions between the giant planets excited eccentricities 
of ~ 0.1, causing an eventual complete depletion in the outer plan- 
etesimal disk and a decreased terrestrial planet formation efficiency 
compared with other stable simulations. Although the planets were 
all roughly one lupiter mass, the two inner planets were actually also 
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Fig. 17. Left: The fraction of systems that would be detectable after 1 Gyr as a function of the eccentricity of the innermost surviving 
giant planet eg for six different wavelengths between 5 and SOO/zm. Each curve represents a horizontal slice through plots such as 
Fig [13] at a different vertical height. The chosen "detection limits" in F IF star were the same as in Fig. [16] 10"** at 5/im, 10"^' at 
15;um, 0.054 at ISjim, 0.55 at IQfim, 10 at 160//m, and 100 at 500/zm. The error bars are 1-cr and were calculated using binomial 
statistics. Note that the different curves are offset for clarity. The bins themselves are logarithmically spaced between 10"^ and 1. 
Right: The fraction of systems with Cg > 0.1 as a function of F/Fs,arO0fim) (1 Gyr). Systems with F/Fsmr < fO"^ are included in 
the bin at F/Fsra,- ~ 10"^. The Spitzer detection limit is shown as the dashed line. This represents a vertical slice through Fig. [T3] 



collisional grinding affects the planetesimals that survive in un- 
stable systems. In some cases the outer planetesimal disk is only 
moderately perturbed by the instability, although in all cases it 
is somewhat disturbed as shown by the lack of unstable systems 
with fluxes as high as the stable systems. After an instability, the 
mass in planetesimals decreases, although the planetesimals' ec- 
centricities and inclinations both tend to increase (not always in 
a simple correlated fashion). Mass loss causes the planetesimal 
population's collisional evolution to slow down and eventually 
stop, as is thought to have occurred in the So lar System's as- 
teroid belt ( Petit et al.ll200Tt iBottke et al.l l2005). Thus, the dust 
flux decreases to a roughly asymptotic value. Once the colli- 
sional timescale for the small particles becomes longer than the 
timescale for Poynting-Robertson our calculation breaks down 
although this occurs at a low enough flux that it should not affect 
our results ( Wvatt et al? '20b7a). 

Figure [18] shows that the debris disk-terrestrial planet corre- 
lation holds in time, especially after 10-100 Myr However, if the 
figure were plotted including the total teiTestrial mass at a given 
time the correlation would hold at even earlier times because in- 
dividual systems would not be restricted to move between panels 
on vertical lines. 



5. Discussion 

In this section we first scale our simulations (section 5.1) to 
match the exoplanet semimajor axis distribution and infer the 
orbital distribution of teiTestrial planets in the cuiTent sample 
(mainly drawn from the radial velocity sample). In section 5.2 
we compare our results with the statistics of known debris disks 
(including cases with known giant planets). In section 5.3 we 



driven slowly across their 2: 1 mutual mean motion resonance after 130 
Myr (which actually decreased the eccentricities). 



discuss the Solar System in the context of our results. In section 
5.4 we discuss the limitations of our simulations. 

5. 1 . Scaling to the observed a-e exoplanet distribution 

The surviving giant planets in our unstable simulations pro- 
vide a match to the observed eccentricity distribution (Fig.|7]l. 
However, in constructing a sample of simulations that represents 
the observed exoplanet systems we do not think that is realistic 
not to include any stable systems for several reasons. First, cur- 
rent attempts to de-bias the observed eccentricity distribution in- 
fer a substantial fraction o f systems - up to ~ 30% - with circu- 
lar or near-circular orbits (IShen & Tumerl2008l ^akamska et al.l 
2010). Second, there exist individual systems that show no ob- 
vious signs of instab ility, for example with planets in reso- 
nance (e.g., GJ 876; iRivera et all 1201 Ol) or wit h many rela- 
tively closely-packed giant planets (e.g., 55 Cancri: lFischer et aLl 
2008). 

Figure [19] (left panel) shows the effect of the contribution 
from stable systems on the goodness of fit to the exoplanet ec- 
centricity distribution. The distribution is best matched when we 
do not include any stable systems in the sample, and drops be- 
low a nominal statistically acceptable limit of p = 0.01 for a 
contribution larger than 17%. We therefore construct our sample 
with a 10% contribution of stable systems, which we think is a 
reasonable compromise between matching the exoplanet eccen- 
tricity distribution and allowing for stable systems. 

The giant planets in our simulations are generally at larger 
orbital distances than the observed giant planets. The observed 
sample of giant planets displays a rapid rise between 0.5-1 AU, 
a plate au to 2 AU, and then a decrease at lar ger orbital dis- 
tances (iButler etal.1 120061: lUdrv & Sant"osll2007h . The rise and 
plateau are real but the decrease beyond 2 AU is an observa- 
tional bias due to the long orbital periods of these planets. The 
actual distribution of giant planets at Jupiter-Satum distances is 
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Fig. 18. The dust-to-stellar flux ratio F/Fsiar at lO/rni as a function of the total mass in terrestrial planets for eight diff'erent times 
from 1 Myr to 3 Gyr The terrestrial planet mass refers to the fi nal value such that simulations move vertically in time on the plot. 
The Spitzer observational limit is shown with the dashed line (iTrilling et al.ll2008h . Grey circles represent stable simulations and 
black circles represent unstable simulations. 



unknown, although observational constraints predict that at least 
10%, and perhaps more than 50%, of stars have giant planets at 
these separations ( Gumming et al. 2008; Gould et al. 2010). 

Our simulations can be scaled to match the combined semi- 
major axis-eccentricity distribution of giant planets beyond 0.5- 
1 AU. At a smaller orbital distance, a giant planet's scattering 
power decreases. This can be quantified by the quantity 0^, the 
ratio of the escape speed from the planet's surfa ce to the escape 
speed from the planet ary system at that location (ISafronovll969t 
iGoldreich et al.ll2004 : 



where Mp and M* are the planetary and stellar mass, respec- 
tively, Rp is the planetary radius and a is the orbital semimajor 
axis. As is inversely proportional to the a, close in giant plan- 
ets scatter less strongly than at larger distances. When 0^ drops 
below 1, collisions become more important than scatteiing. For 
our simulations, only the lowest-mass giant planets would drop 
to 0^ < 1 by scaling them to 0.5-1 AU. In addition, the planet- 
planet scattering mechanism has been proven at this range of 



distances (iMarzari & WeidenschillingI 120021 lJuric & Tremaind 
2008). Thus, we conclude that it is dynamically appropriate to 
scale the giant planets inward to 0.5-1 AU. We cannot, however, 
scale to closer distances because the dynamical regime is diff'er- 
ent and giant planet-planet collisions are more likely than scat- 
tering events. 

We assume that the underlying distribution that is being 
probed by current (mainly radial velocity) observations rises lin- 
early from zero at 0.5 AU to 1 AU, where it flattens off and 
remains constant to 5 AU. We draw randomly from this distri- 
bution and scale the innermost giant planets from ten randomly 
chosen simulations - nine that were unstable and one that was 
stable according to relative contribution of stable vs. unstable 
systems in our sample - to match this value. We then re-scale 
the teiTestrial planets in each system to the same size scale as 
the inner giant planet's orbit. Our sample was created by repeat- 
ing this 100,000 times with different random numbers. 

Figure |20] shows the radial distiibution of terrestrial and gi- 
ant planets after this scaling. In essence, this figure shows the 
expected radial distribution of terrestrial planets in the known 
extra-solar giant planet systems. In these systems, our simula- 
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Fig. 19. The p value from a Kolmogorov-Smirnov test compar- 
ing the eccentricities of our simulations vs. the observed exo- 
planet distribution as a function of the fraction of stable systems 
included in the sample. The p value drops below an acceptable 
level of 0.01 for more than a 17% contribution from stable sys- 
tems. 
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Fig. 20. Semimajor axis distribution of simulated terrestrial plan- 
ets (dashed line) from our set of simulations, derived by scaling 
the innermost surviving giant planet in each simulation to match 
an assumed underlying distribution for relevant exoplanets that 
increases linearly from zero at 0.5 AU and is constant from 1-5 
AU. 



tions predict a factor of a few higher abundance of terrestrial 
planets at a few tenths of an AU than at 1 AU because, given the 
typical giant-terrestrial planet spacing, the formation of planets 
at 1 AU requkes distant giant planets that are hard to detect by 
current methods. Planets within ~ 0.1 AU are sparsely populated 
because of the assumed inner edge of the embryo disk at 0.5 AU. 
The peak in the frequency of terrestrial planets at a few tenths of 
an AU is due to a combination of our inner disk edge and the 
giant planets' radial distribution. 

As noted above, the scaling we performed is dynamically 
permissible: it does not change the regime of accretion of the 
terrestrial planets nor the scattering regime of the giant planets. 
However, our simulations each contained a constant mass in em- 
bryos and planetesimals in the inner disk. By scaling the giant 



and terrestrial planets, we are effectively re-scaling the initial 
disk masses such that the same amount of mass would initially 
have been placed into an annulus whose position and width can 
vary. In addition, if we had included initially closer-in material 
in our simulations, the peak in Fig|20]would have shifted inward. 
Despite these limitations, in the regime that we have considered 
we think that the shape of the curve is physical, and we predict 
that, at least within the known sample of extra-solar giant plan- 
ets, terrestrial planets at ~ 0.3 AU should be several times more 
abundant than at 1 AU. 



5.2. Comparison between our simulations and observed 
debris disks 

Observations (mainly with NASA's Spitzer Space Telescope) 
have shown that roughly 15% of solar-type stars younger 
than 3 00 million years have measurable dust fluxes at 
24pm (Gaspa r et al.l l2009h but that this fraction decreases in 
time and flattens off at -3% JCarpenter e t al. 2009). At 70/im, 
16% of stars observable du st and there is no m e asured decrease 
in this fraction with age (iTrilUng et al.l 120081: [Carpenter et al.l 
2009). Considering the current exoplanet/debris disk sample, 
9% + 3% of planet-hosting stars were detected at 24 or 70A<m 
comp ared with 14% + 3% for stars without planets (Brvde n et al.l 
|2009|) . An update of that study found debris disks around the 
same fracti on, ~ 15%, of star s with and without known gi- 
ant planets dKospal et alj 120091) . In addition, the strong corre- 
lation be tween stellar metallicity and the fraction of stars with 
plane ts (iGonzalezI [T997I: ISantos et al.ll200U iFischer & Valentil 
l2005h does not hold for the current sample of debris disks, 
whos e presence is metallicity-independent ("Beichman et al 



20061: iGreaves et al.l 12006: Brydenetal. 2006; Kospal et aL 



20091) . 



Figure[T3]shows that, in our simulations, the dust flux is anti- 
correlated with the giant planet eccentricity. Almost all lower- 
eccentricity (e < 0.1 - 0.2) giant planets are in systems with 
debris disks, but at higher eccentricities the fraction of dusty 
systems decreases as does the dust brightness itself. Figure [17] 
shows that the fraction of systems that are detectable at all wave- 
lengths from 5pm to 500/im decreases with increasing eccen- 
tricity of the innermost surviving giant planet. In addition, the 
fraction of planets with eg > 0.1 decreases with F/FsiarUOl-iin), 
including a dramatic drop for F IF star > 10. 

We therefore expect to see a correlation between the orbital 
properties of exoplanets and th e detectability of cold dust. The 
dataset of iBrvden et al.l (l2009l) detected debris disks at 70;um 
around 13 planet-hosting stars out of 146 for a detection fre- 
quency of 8.9%;^j g^°. We arbitrarily divide the sample in two 
based on the eccentricity of the innermost giant planet in each 
system at a value of 0.2. Debris disks are detected in 8 of 76 
systems (10.5%;^26%)i'^ '■^^ low-eccentricity subsample and 5 
of 70 systems (7.1%^'f g*) in the high eccentricity subsample. 
Although the detection rate is somewhat higher in the lower- 
eccentricity subsample, the difference is not statistically signifi- 
cant and we must wait for better statistics from larger surveys. 

Our simulations do produce some systems with high- 
eccentricity giant planets and bright dust emission (Fig. [T3T l. in 
agreeme nt with the detected deb ris disks in known exoplanet 
systems (iMoro-Martm et an2010h . In these cases the dynamical 
instability tends to be asymmetric and confined to the inner plan- 
etary system and these are therefore not generally good candi- 
dates for terrestrial planets, which also agrees with the observed 
systems (IMoro-Martm et afl2010l) . The outcome of a given sys- 
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tern depends critically on the details o f the instability, which is 
determined by the giant planet masses (iRavmond et alJlToiOl) . 

Our approach is not unique; other approaches based on dust 
production by collision al cascades can also reproduce the de- 
bris di sk observations ('Krivov et al.' '2005', ' 20061: IWyatt et al 



2007bl: FW vatt 2008; Kenvon & Bromley 200: 



uua IWy 
181 120101: 



Krivov 



2O10H Kennedv & Wvatt 2010). In addition, there is a qual 



itative difference between models that consider planetesimal 
disks to be "self-stirred" (i.e., eccentricitie s are excited by ac- 
creting bodie s within planetesimal disks; 'Kenvo n & Bromlevl 
[2008, 2010; Kennedv & Wyatt 2010) or those that consider ex - 
ternal sources fo r planetesimal stiiTing (e.g. Wvatt et al. 2010*). 
iMustill & WvattI (1200 9) showed that self-stiiTed planetesimal 
disks tend to be fainter than disks stirred by giant planets. Debris 
disks may be explained by some combination of these ideas. 

Nonetheless, we conclude that our simulations are consis- 
tent with the known sample of debris disks in exoplanet systems. 
However, our initial conditions are biased in that all systems that 
could potentially produce debris disks also contain giant planets, 
which is not consistent with the observation that the frequency 
of debris disk + giant exoplanet systems is a bout the same as de- 
bris disks with no d etected giant exoplanets jBrvden et al.ll2009t 
iKosnal et alJl2009h . In paper 2 we use multiple sets of simula- 
tions to construct a sample that adequately matches the entire 
debris disk and exoplanet samples, and use that sample to infer 
the properties of as-yet-undetected terrestrial exoplanets. 

5.3. Our Solar System in context 

Our results suggest that the Solar System is unusual at the ~ 
15-25% level. This corresponds to the fraction of simulations 
that form three or four terrestrial planets (Fig. including a 10% 
contribution from stable systems. By the same weighting 38% of 
simulations destroy all their terrestrial material.). 

The Solar System lies at the very edge of the debris disk 
correlations in Figs. [T3] and [T4]because of its combination of a 
rich terrestrial planet system, a low-eccentricity innermost gi- 
ant planet, and a low dust flux. To a distant observer, the Solar 
System's faint debris disk would suggest a dust-clearing insta- 
bility in the system's past. However, Jupiter's low-eccentricity 
orbit would imply that the instability was weak and that the sys- 
tem may in fact be suitable for terrestrial planets J3 This naive ar- 
gument is remarkably consistent with our current picture of the 
LHB instability as an asymmetric, outward-directed instability 
that included a scattering event bet ween Jupiter and an ice giant 
but not between Jupiter and Saturn (Morbidelli et al.ll2010l) . 

The Earth provides an interesting test case. On long 
timescales Earth's eccentricity oscillates between 0.0002 and 
0.058 and its inclination between zero and 4.3° dOuinn et al.l 
[T99TllLaskar" et al.|[T993l) . When compared with the stable sim- 
ulations, the Earth's time-averaged eccentricity is significantly 
smaller than the median and its inclination is also smaller 
However, Earth's oscillation amplitudes are more than 50% 
larger than the median values for the stable systems. This sit- 
uation is the same for Venus' orbit. This presents a confusing 
situation; perhaps Earth and Venus' e and / are lower than these 
simulations because we are limited in the numerical resolution 
needed to adequately model dissipative processes. But if that 
were the case, we would expect Earth and Venus' oscillation 
amplitudes to also be smaller than the simulated planets'. One 
explanation for the Earth and Venus' orbits is that they formed 



in a dissipative environment but that were later dynamically per- 
turbed during the instability that caused the late heavy bombard- 
ment (Gomes etal. 2005; Brasseretal. 2009; Morbidelli et aT| 
12010 ). The perturbation was not large enough to disrupt the sys- 
tem's stability but sufficient to increase the amplitude of orbital 
oscillations of Earth and Venus. 

The instability that caused the LHB is not captured in our 
simulations nor in current exoplanet observations. The degree 
to which our simulations interpret that the Solar System is un- 
usual depends on how well the simulations characterize the in- 
stabilities in extra-solar planetary systems. As our simulations 
reproduce the observed exoplanet eccentricity distribution with 
no free parameters, it appears that our simulations do in fact 
capture the essence of instabilities among the known exoplan- 
ets. However, an instability such as the one proposed by the 
Nice model leaves little to no trace because the giant planets' 
eccentricities remain very small (ej and es are only ~ 0.05). 
It is plausible that Nice model-type instabilities are common in 
outer planetary systems, although if they systematically destroy 
their outer planetesimal disks then debris disk statistics con- 
strains the fraction of stars that undergo such instabilities more 
than about 10 Myr afte r stellar formation to be less than about 
10% (iBooth et al.ll2009l) . In Paper 2 we present an additional set 
of simulations with larger mass ratios between the giant planets 
in which Nice model-type instabilities can occur Such instabili- 
ties do not change our conclusions; Nice model instabilities can 
effectively be lumped in with the stable systems as their impact 
on inner planetary systems is smallj^ 



5.4. Limitations of our approach 

As with any numerical study, our simulations are limited in sev- 
eral ways. 

Our most important assumption is that the inner and outer 
regions of protoplanetary disks are connected such that obser- 
vations of debris disks in the outer parts of planetary systems 
can tell us something about terrestrial planets in the inner parts 
of these systems. However, there exist substantial uncertainties. 
For example, several relevant processes - notably the formation 
of planetesimals as well as giant planets - are only modestly- 
well understood. In addition, it is unknown whether the efficien- 
cies and timescales of these processes vary with distance from 
the star If there exists a systematic bias to create an imbalance 
between the inner and outer disk mass, it could qualitatively af- 
fect our results. For example, if planetesimal formation is much 
more efficient in the outskirts of planetary systems then the typ- 
ical system may contain an outer planetesimal belt but no inner 
planetesimals or embryos from which to form tetTestrial planets. 
Alternately, one can imagine that outer planetesimal disks might 
be systematically destroyed in systems that do form terrestrial 
planets, as was the case for the Solar System. If one of these 
scenarios is true, then the initial conditions for inner and outer 
planetary systems may not be coupled as strongly as we have 
assumed, and their outcomes may not cotTelate as strongly as in 
our simulations. 

Despite these uncertainties, we think that our approach, 
in particular the assumption that inner and outer disks are 
connected, is the simplest interpretation of cutTent observa- 



Of course, it would take about a decade of precise observations for 
these aliens to pin down Jupiter's orbital eccentricity. 



' iBrasser et al ] (l2009h showed that some Nice model instabilities can 
destabilize the orbits of the terrestrial planets due to sweeping secular 
resonances and potentially cause collisions between the terrestrial plan- 
ets. However, this process does not remove all the terrestrial planets and 
so, in the context of our results, nothing has changed. 
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tions and theory. The frequency of close-in planets is 12% 
in the 3-10 range and, by extrapolation, 23% in the 0.5- 
2 Me range (Howard etal. 2010). The ob served freque ncy of 
debris disks around FGK stars of 16% dTrilling et all l2008h 
represents a lower limit for the frequency of outer planetesi- 
mal disks. Preliminary results from the Herschel DUNES sur- 
vey ( Eiroa et al. 2010) suggest that roughly 1/3 of stars have de- 
bris disks (C. Eiroa, personal communication), which means that 
that the frequency of inner planets and outer planetesimal belts 
are probably within a factor of a few or less. Although this cer- 
tainly does not prove that our initial conditions are correct, it 
does provide circumstantial support for our basic assumption of 
a connection between the inner and outer disk although we note 
that there is as yet no observed correlation between the two. 

Our initial conditions, though chosen to match models of ear- 
lier phases of planetary growth, are ad hoc. All planetary systems 
in our sample contain the same mass in terrestrial embryos and 
planetesimals (9Me), form their innermost giant planet at 5.2 
AU, and contain the same mass in an outer disk of leftover icy 
planetesimals (50 M®). In reality, there is a spread of several or- 
ders of magnitude in the disk mass (e.g., Andrews & Williams 
l200 7a) that affects both the types of pl anets that form and 
the amount of leftove r planetesimals (e.g.. iGreaves et al" l l2007t 
iThommes et"al]|2008h . In addition, rec ent observations sugg est 
that low-mass disks are more compact ('Andrews et al."20ldl) so 
there may be a correlation between the dis k mass and the loca- 
tion of planet formation as well (see also iKennedv & Kenvoiil 
I2OO81) . The disks that we modeled are comparable in mass 
to the minimum-mass solar nebula model and are probably 
more mas sive than the typical disk dEisner & Carpenter! l2003i 
lAndrews & Williams 2007a). 

Our simulations are confine d to Solar-mass stars, w hich are 
a small minority of all stars (i Bochanski et al.l l2010h . To ex- 
pand on this work we should consider additional stellar types. 
Debris dis ks are c urrently very difficult to detect around low- 
mass stars (iGautier et al. 2007) but are extremely interesting as 
planet hosts because they dominate the stellar population of the 
Galaxy and are very long-lived. In contrast, debris disks are 
much easier to detect around A stars, but these are relatively 
few in number and their lifetimes are much shorter There may 
be interesting differences in the evolution of planetary systems 
around other stellar types. 

There are several physical processes not included in our 
simulations. In particular, we did not include the effects of gi- 
ant planet migration because population synthesis models are 
currently unable to repr oduce the observed exoplanet mass 
and orbital distribution (iHoward et al.l l2010l) . Including mi- 
gration would have the benefit of providing a natural trig- 
ger for giant planet ins tabilities (Adams & Laughlin 2003] 
iMoorhead & Adamsll2005b . althou gh this depends on the details 
of the depletion of the gaseous disk d Crida et al.l2008h . However, 
given that the giant planet observations can be matched with lit- 
tle to no giant planet migration, we chose not to include it. 

6. Conclusions 

Our main results are as follows: 

- Giant planet instabilities are destructive to terrestrial 
planet formation. The survival of terrestrial embryos and 
planetesimals depends on the strength of the instability 
as measured by the minimum giant planet perihelion dis- 
tance (Fig[8j see also lVeras & Armitag92006l) . In about 40% 
of our unstable simulations all rocky material was removed 



from the system, in large part by being thrown into the host 
star. About 1/5 of our unstable simulations produced a sys- 
tem containing a single terrestrial planet (Fig.|9). 

- Terrestrial exoplanets on excited orbits should be com- 
mon. The median eccentricity of surviving terrestrial plan- 
ets in our simulations was about 0.1, but the distribution 
extends above 0.5 (Fig. |7|. The most excited orbits belong 
to single-terrestrial planet systems. Compared with systems 
with many terrestrial planets, single-planet systems have 
only slightly higher eccentricities and inclinations but their 
oscillations in e and / are far larger (Fig. [TOi . 

- Debris disks are anti-correlated with strongly-scattered 
giant planets. Strong scattering events produce eccentric gi- 
ant planets with large radial excursions that dynamically de- 
plete the outer planetesimal disk by exciting their orbits until 
they cross a giant planet's at which point they are quickly 
ejected from the system. Thus, we expect continued obser- 
vations to show an anti-correlation between the fraction of 
systems with debris disks and the giant planet eccentricity. 

- Debris disks correlate with a high efficiency of terrestrial 
planet formation. Strong scattering events yield large gi- 
ant planet eccentricities, and these eccentric giant planets 
tend to disturb both the inner and outer planetary system. 
Thus, in a strongly perturbed system the giant planets tend 
to destroy both terrestrial planetary embryos - aborting the 
growth of terrestrial planets - and outer planetesimals - pre- 
venting the creation of debris disks by long-term collisional 
evolution. In contrast, in a calm system the giant planets will 
not strongly impede on the inner or outer planetary system, 
allowing for the formation of terrestrial planets and long- 
lasting cold dust. The debris disk-terrestrial planet correla- 
tion holds for all wavelengths we tested but is clearer for 
A > 15/im (Figs. [TSlandfTSIl. The correlation also holds for 
all times later than about 10-30 Myr (Fig. [18) and probably 
even earlier 

- Within the known sample of extra-solar giant planets, 
terrestrial planets at a few tenths of an AU should be sev- 
eral times more abundant than either terrestrial or giant 
planets at 1 AU. In section 5.1 we scaled the outcomes of 
our simulations to match the observed semimajor axis and 
eccentricity distributions of giant exoplanets. This scaling 
produced a radial distribution of terrestrial exoplanets that 
peaks at a few tenths of an AU and drops below the giant 
planet frequency at 1.3 AU (Fig.l20t. However, we note that 
the distribution is incomplete in its inner regions due to our 
initial conditions, in particular the inner edge in our embryo 
distribution at 0.5 AU. 

- The Solar System lies at the outskirts of several correla- 
tions, probably because of the instability that caused the 
Late Heavy Bombardment. The Solar System has a rare 
combination of multiple terrestrial planets, low-eccentricity 
giant pla nets, and very low dust content (currently, F/Fstar ~ 
2x 10~^: lBooth et alj20"09 '). In addition. Earth and Venus' or- 
bits are more circular and coplanar than terrestrial planets in 
typical stable planetary systems but they have significantly 
larger-amplitude oscillations in those quantities. This can be 
explained if the Solar System's formation was quiescent but 
it underwent a later punctual event that was strong enough 
to remove most of the outer planetesimal disk and give the 
inner Solar System a small kick but did not destabilize the 
inner Solar System or impart a large eccentricity to Jupiter. 
This is in agreement with the current picture of the instabil- 
ity that caused the LHB (Morbidelli et al. 2010). This type 
of instability is poorly constrained by observations, is much 
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weaker than the instabilities inferred from the exoplanet ec- 
centricity distribution, and is not captured in the simulations 
presented here. 

One implication of our results is that exo-zodiacal dust 
clouds around stars with terrestrial planets may often be brighter 
than the Solar System's zodiacal cloud. A major obstacle to 
the direct imaging of terrestrial exoplanets is the amount of 
bright d ust close to t hose planets, i.e., their exo-zodiacal dust 
clouds dCashl l2006t iDefrere et al.1 120101: iNoecker & Kuchnerf 
I2OIOI) . The Solar System's zodiacal dust has been shown to de- 
rive almost entirely from Kuiper Belt comets that were scattered 
by the giant planets into the inner Solar System, where they par- 
tially sublim ated to produce w arm dust before eventually being 
ejected ( Nesvomv et al.llToiOl) . Around other stars, cold debris 
disks should trace the same population of comets that produces 
exo-zodiacal dust: debris disks represent planetesimals on sta- 
ble orbits in the outer system and exo-zodiacal dust is gener- 
ated by the small fraction of bodies that has been destabilized 
and is in the process of being removed from the system. Our 
results suggest that systems with bright debris disks are excel- 
lent targets in the search for terrestrial exoplanets. These sys- 
tems contain at least a few (and often more than 20Me) in 
surviving cometary mat erial, 1-2 orders of ma gnitude more than 
the current Kuiper Belt jGladman et al.ll200lT) . If the comet flux 
scales with the number of outer planetesimals then systems with 
bright debris disks should also harbor bright exo-zodiacal clouds 
close to the terrestrial planet zone. However, the dynamics of 
the outer planetary systems - in particular the architecture and 
masses of the giant planets - are key in determining the fluxes of 
new comets in these systems as well as their residence lifetimes 
in the inner planetary systems. In addition, there is almost cer- 
tainly a significant population of systems with terrestrial plan- 
ets without bright debris disks, i.e. what we have called "false 
negatives". Those systems may be harder to find because they 
lack debris disks to signpost the presence of terrestrial planets, 
but they could prove easier to image because they may contain 
much fainter zodiacal clouds. We plan to test these ideas in fu- 
ture work. 

In a companion paper (Raymond et al 201 1 ; referred to in the 
text as paper 2) we explore the effect of several other parameters 
on the results obtained here. In particular, we present results of 
several other sets of simulations that vary the giant planets' mass 
distribution and total masses, the width of the outer planetesimal 
disk, the existence of icy embryos within the outer planetesimal 
disk, and the presence of disk gas at the time of giant planet 
instabihties. In that paper we confirm the large-scale results pre- 
sented here but with several important clarifications and depen- 
dencies. We also carefully match our simulations to the observed 
statistics of giant exoplanets and debris disks to obtain an esti- 
mate for the fraction of stars that host terrestrial planets, rjEarth 
in the famous Drake equation. 

In a second companion paper (Raymond et al 201 lb), we ex- 
plore the fate of bodies - planetesimals, planetary embryos, and 
giant planets - that are ejected from unstable planetary systems 
and pollute the galaxy. By matching giant exoplanet and debris 
disk statistics we estimate the abundance of this population of 
free-floating bodies and the chances for their unambiguous de- 
tection either in interstellar space or entering the Solar System. 

7. Appendix A: Numerical Tests 

We performed simple numerical tests to choose an inner bound- 
ary appropriate for this timestep by placing a particle at 1 AU on 




q(AU) 

Fig. 21. Fractional error in the semimajor axis a of a Mars-mass 
particle with initial a - \ AU as a function of perihelion distance 
q. The particle was initially placed on a circular, nearly polar 
orbit (initial inclination of 89.9°) in the presence of an outer giant 
planet. The Kozai effect forced the particle to fall into the star, 
and we tracked the integration error in time. 

an orbit that was highly inclined with respect to a distant, Jupiter- 
mass plane{3- In this configuration, the Kozai effect (iKozail 
1962) drives a dramatic increase in the particle's eccentricity. By 
tracking the integrator error at ever-smaller perihelion distance, 
q - a{\ — e), we saw that the fractional error in the semimajor 
axis dala increased to greater than lO"** inside roughly 0.2 AU, 
as shown in Figure |2T| We therefore chose 0.2 AU as our inner 
particle boundary, inside which objects are removed from the 
simulation and assumed to have collided with the star We chose 
100 AU as our corresponding outer boundary, beyond which 
bodies are assumed to have been dynamically ejected from the 
system. 

With our choice of timestep and inner boundary, the majority 
of our simulations maintained excellent (dE/E < 10^*) energy 
and angular momentum conservation properties^. However, as 
shown in Figure |22] a fraction of the simulations in which the 
minimum giant planet perihelion distance was small {q <2 AU) 
exhibited substantially worse conservation properties. Given this 
behavior, we adopted an empirical energy conservation thresh- 
old of dE/E < 1 X 10"^, and discarded from the final sample any 
runs that failed to meet this limit. 

Does our cutoff in energy at dE/E < 1 x 10"^ bias our re- 
sults? The true outcome of our numerical experiment depends on 
the details of the giant planets' orbital evolution, and it is only 
relatively extreme cases with q < 2 AU for which significant in- 
tegration error occurs. It is for these close perihelion passages 
that all terrestrial material is destroyed, called mode 2 accre- 
tion in the main text. In the case of the five simulations with 
dE/E > 1 X 10"^, two contained a single surviving terrestrial 
planet and three had destroyed all of their terrestrial planets. By 
removing these cases we are therefore slightly biasing our sam- 
ple at the 5% level away from accretion modes 2 and 3, and so 
we include this small extra contribution in our estimates in the 
paper of the fraction of systems for which the different accretion 
modes occur 

We thank Hal Levison for suggesting this numerical test. 
" We are aware that removing particles from the simulation, at any 
radius larger than the physical one (approximately the size of the star), 
can in principle result in unphysical behavior, even when energy and 
angular momentum are well-conserved. Unfortunately, it is currently 
infeasible to run long-duration terrestrial planet formation simulations 
with dramatically shorter timesteps and smaller inner boundary radii. 
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Fig. 23. Eccentricity of the innermost giant planet (left panel) and dust-to-stellar flux ratio at 70yum after 1 Gyr of evolution (right 
panel) as a function of the fractional error in the system's energy budget dE/E. This plot is only for systems with minimum giant 
planet perihelia of less than 2 AU (see Fig. l22l i. Our energy error cutoff at 0.01 is shown with the dashed line. There is no evidence 
for contamination of our sample. 
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Fig. 22. Fractional error in the system's energy budget dE/E as 
a function of the smallest perihelion distance of a giant planet 
for all simulations. 

We do not see any clear signature of other bias introduced 
in our analysis from our energy cutoff. To assess this possibility, 
we restrict ourselves to systems for which the minimum giant 
planet perihelion was less than 2 AU because this is where large 
errors occur. Figure 22] shows the giant planet eccentricity and 
the dust-to-stellar flux ratio at 70/im after 1 Gyr as a function 
of dE/E for this subsample. Both panels are scatter plots, with 
no clear trend or any indication that the computed dE/E value 
changes the outcome in any way. We therefore conclude that our 
chosen cutoff in integrator energy, while less stringent than some 
other studies, has no measurable impact on our results. 
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